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Preface 


The  purpose  of  this  study  was  to  explore  the  possible 
application  of  a  thrust  vector  control  law  constraining 
either  the  radius  of  apogee  or  perigee  to  transfers  between 
circular  and  eccentric  orbits  for  continuous,  lew  thrust 
spacecraft.  I  have  tried  to  include  as  many  details  as 
possible  in  the  analysis  to  aid  anyone  who  may  want  to 
perform  similar  or  follow-on  research. 

To  prove  the  validity  of  the  assumptions  made  in  this 
study,  I  thought  it  necessary  to  perform  a  separate 
performance  analysis  of  the  more  promising  electrical 
propulsion  systems  (Appendix  A) .  I  am  truly  grateful  to  Mr 
Mike  Patterson  (NASA  Lewis) ,  Capt  Wayne  Schmidt  (AFAL) ,  and 
fellow  Boilermaker,  Mr  Joe  Cassady  (Rocket  Research  Company) 
for  providing  me  with  the  data  I  needed  to  put  together  a 
reasonably  valid  performance  estimate  of  an  Ammonia  arejet 
and  Xenon  ion  propulsion  system. 

My  sincere  thanks  goes  to  Dr  William  Wiesel,  my  thesis 
advisor,  for  his  assistance  and  patience  over  the  past 
several  months.  Also,  I  thank  my  wife,  Mary,  and  my 
children,  Bethany,  Dustin,  and  Brandon,  for  their  under¬ 
standing  and  support  during  those  many  hours  that  I  spent 
away  from  them  over  the  past  year. 

I  dedicate  this  thesis  in  memory  of  my  grandparents,  the 


late  Ora  and  Rita  Beeker. 
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Abstract 

The  purpose  of  this  study  was  to  investigate  the  Av 
requirements  of  a  continuous,  low  thrust  spacecraft  perform¬ 
ing  coplanar  orbit  transfers  that  are  constrained  by  either 
constant  radius  of  apogee,  perigee,  or  a  combination  of  the 
two.  The  transfers  were  separated  into  two  timescale 
problems.  The  fast  timescale  involved  an  optimization  of 
the  planar  thrust  control  angle,  a,  to  produce  the  maximum 
change  in  eccentricity  or  semimajor  axis  over  a  single 
revolution.  The  slow  timescale  applied  the  fast  timescale 
results  to  complete  the  transfer  through  many  revolutions 
about  the  primary  body.  The  constrained  radii  control  laws 
developed  provide  optimal  circular-to-eccentric  and 
eccentric-t o-eccentric  orbit  transfers.  However,  when 
applied  to  circular-t o-circular  transfers,  the  resulting  Av 
is  nearly  twice  that  obtained  using  the  most  optimal 
continuous  thrust  control  law  (a  =  0),  i . e .  "spiral" . 

Future  recommended  studies  include  the  development  of 
control  lavs  to  provide  specified  changes  in  apogee, 
perigee,  and  the  argument  of  periapsis. 
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CONTINUOUS,  LOW  THRUST  COPLANAR  ORBIT  TRANSFERS 
WITH  VARYING  ECCENTRICITY 


1_.  Introduction 

The  United  States'  decision  to  construct  a  space 
station  within  the  next  decade  has  launched  this  country 
into  a  new  era  in  which  it  will  be  accessing  space  more  than 
ever  before.  With  this  new  era  comes  an  increasing  need  for 
an  energy  efficient  orbit  transfer  vehicle  to  act  as  a 
freighter  between  low  and  high  earth  orbit  -  a  vehicle  vital 
in  reducing  the  work  loads  of  supporting  expendable  boosters 
and  the  Space  Transportation  System. 

A  transfer  vehicle  equipped  with  a  continuous,  low 
thrust,  electric  propulsion  system,  such  as  the  ion  systems 
currently  in  development  at  NASA  Lewis  Research  Center  or 
the  arcjet  systems  at  the  Air  Force  Astronautics  Laboratory, 
could  provide  for  this  need,  greatly  reducing  fuel  and 
support  requirements  (4;  6;  10:457-467) .  However,  low 
thrust  systems  of  this  type  require  many  revolutions  to 
complete  their  transfer  and  are  therefore  limited  to 
missions  where  extended  transfer  times  (months)  are  accept¬ 
able.  In  addition,  the  electronics  of  the  transfer  vehicle 
and  payload  must  be  able  to  endure  the  extended  exposure  to 
the  high  radiation  within  the  lower  Van  Allen  belt. 


Alfano  Cl)  addressed  the  problem  of  locating  an  optimal 
thrust  profile  for  noncoplanar  transfers  performed  by  a 
continuous,  low  thrust  vehicle  between  circular  orbits. 
However,  the  more  complex  problem  of  locating  an  optimal 
thrust  profile  for  transfers  between  eccentric  orbits  has 
never  been  addressed,  although  there  are  many  useful 
applications,  such  as,  the  placement  of  communications 
satellites  into  Molniya  orbits. 

To  provide  an  initial  investigation  into  the  eccentric 
transfer  problem,  this  research  examines  coplanar,  eccentric 
transfers  constrained  by  holding  the  radius  of  apogee  or 
perigee  constant  over  each  orbit.  The  development  of  this 
problem  is  simplified  by  dividing  the  derivation  between  a 
fast  and  slow  timescale,  similar  to  that  done  by  Alfano. 

The  fast  timescale  problem  involves  an  optimization  of  the 
change  in  the  orbital  elements  over  one  revolution  while 
implementing  the  constraint  relation  and  holding  the  vehicle 
mass  and  acceleration  constant.  The  slow  timescale  problem 
combines  the  fast  timescale  results  with  the  changing  mass 
and  acceleration  to  complete  the  transfer  over  many 
revolutions . 

The  development  of  the  fast  timescale  solution  in 
chapter  two  shows  this  constraint  relationship  is  unique 
since  it  provides  two  control  laws  for  the  planar  thrust 
angle,  a,  (one  for  maximizing  Aa  and  the  other  for 
maximizing  Ae)  which  result  in  identical  changes  to  the 


orbit’s  semimajor  axis  and  eccentricity.  Thus,  either  of 
these  two  control  laws  are  applicable  to  any  coplanar 
transfer,  whether  the  transfer  involves  a  change  in 
eccentricity,  the  semimajor  axis,  or  both.  However,  these 
results  Cand  the  slow  timescale  problem  of  chapter  three) 
indicate  that  utilizing  these  control  laws  for  transfers 
between  coplanar  circular  orbits  is  much  less  efficient  than 
the  spiral  (a  =  0)  control  law. 
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II .  The  Fast  Timescale  Problem 


Pr ob 1 em  Stat  ement 

The  fast  timescale  analysis  performed  by  Alfano  (1:3-23) 
included  an  optimization  of  the  out-of-plane  thrust  angle, 

9,  for  each  of  the  two  orbital  parameters  involved,  the 
change  in  semimajor  axis,  Aa ,  and  the  change  in  inclination, 
Ai ,  maximizing  one  parameter  for  any  given  value  of  the 
other.  Optimization  problems  involving  transfers  between 
eccentric  orbits  become  more  complex  due  to  the  addition  of 
the  planar  thrust  angle,  a,  eccentricity,  e,  and  argument  of 
periapsis,  w.  Thus,  the  fast  timescale  analysis  includes  an 
optimization  of  both  a  and  9  for  each  of  the  four  parameters 
involved,  Aa ,  Ai ,  Ae,  and  Aw. 

This  problem  can  be  simplified  by  considering  only 
coplanar  transfers  and  by  placing  no  restrictions  on  Aw. 

This  results  in  an  optimization  of  a  only  (since  0=0)  for 
the  parameters  Aa  or  Ae  (Ai  =  0) .  However,  the  application 
of  the  resulting  control  law  would  produce  uncontrollable 
changes  in  the  radius  of  perigee,  possibly  allowing  an 
impact  with  the  primary  body.  Therefore,  the  problem 
constraints  should  include  the  radius  of  perigee  instead  of 
the  eccentricity  or  semimajor  axis. 

The  coplanar  transfer  problem  addressed  by  this  study 
is  constrained  by  holding  the  radius  of  apogee  or  perigee 
constant.  Thus,  the  optimization  is  further  simplified 


since  the  constraint  relationship  is  no  longer  dependent  on 
the  total  transfer  and  allows  the  Lagrange  multiplier,  X,  to 
be  found  in  the  fast,  rather  that  the  slow  timescale 
problem.  Given  these  constraints,  the  fast  timescale 
problem  involves  an  optimization  of  a  in  order  to  maximize 
the  magnitude  of  Aa  or  Ae  for  one  orbit  about  the  primary 
body  (Earth) . 

The  formulation  of  the  fast  timescale  problem 
is  based  on  the  following  three  assumptions: 

First,  the  problem  is  restricted  to  two  bodies, 
the  earth  and  the  spacecraft,  which  are  modeled  as  point 
masses . 

Second,  since  electric  propulsion  systems  have 
very  low  propellant  mass  flow  rates,  the  percent  change  in 
mass  of  the  spacecraft  over  one  orbit  is  small.  Thus,  the 
spacecraft's  mass  will  be  considered  constant.  Appendix  A 
provides  data  in  support  of  this  assumption,  showing  that  a 
5000  kg  transfer  vehicle  system  (vehicle  and  payload)  has  a 
change  in  mass  over  one  orbital  period  of  less  than  3%  if 
equipped  with  the  proposed  arcjet  propulsion  system  and  less 
than  .25%  if  equipped  with  the  ion  system. 

Finally,  the  low  thrust  associated  with  electric 
propulsion  systems  produces  such  small  accelerations  on  the 
spacecraft  that  the  changes  in  the  orbital  elements  occur 
very  slowly.  This  leads  to  the  assumption  that  these  these 
parameters  vary  so  little  during  the  period  of  one  orbit 
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that  they  can  be  treated  as  constants.  Appendix  A  also 
provides  justification  by  showing  the  maximum  acceleration 

—3  2 

as  seen  by  the  proposed  TVS  is  approximately  4.5  x  10  m/s 

-*-3  2 

using  the  arcjet  system  and  1.3  x  10  m/s  using  the  ion 
system . 

Derivation 

Coordinat e  System  Definition.  The  coordinate  system 
utilized  throughout  this  analysis  (2:397,398)  is  located  at 
the  center  of  mass  of  the  transfer  vehicle,  with  its 
principle  axis,  R  (unit  vector  R)  ,  located  along  the 
instantaneous  radial  vector,  r.  The  second  axis,  S  (unit 
vector  S) ,  is  located  in  the  orbital  plane  perpendicular 
to  R  in  the  direction  of  increasing  true  anomaly,  v .  The 
third  axis,  W  (unit  vector  W) ,  is  perpendicular  to  both  R 

AAA 

and  S,  forming  an  orthogonal  triad  where  W  s  R  x  S.  Thus, 
this  coordinate  system  is  rotated  by  an  angle  v  with  respect 
to  the  perifocal  coordinate  system.  Figure  2-1  defines  the 


geometric 
perifocal 
(2  : 5  3-5  9) 


relationships  among  the  spacecraft  (RSW) , 
(PQW) ,  and  geocentric-equatorial  (IJK)  frames 
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direction  1 


Teroal  eqinoz 
direction 


Figure  2-1 . 


Relationships  among  the  Spacecraft , 
Perifocal,  and  Geocentric-Equatorial 
Coordinate  Frames 
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Vehicle  Thrust  Vector .  The  thrust  vector,  T,  as 
defined  in  the  RSW  coordinate  system,  is  directed  at  an 
angle  9  with  respect  to  the  orbital  plane  formed  by  the 
radial  unit  vector,  R,  and  the  tangential  unit  vector,  S 
(Figure  2-2) .  The  projection  of  T  on  the  orbital  plane 
is  located  at  an  angle  a  with  respect  tc  S.  The  acceleration 
vector  lies  along  the  thrust  vector  and  is  defined  as 


T 

m 

TVS 

A 

T 

<« 

b 

« 

III 

+ 
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0  <  a  <  2rc;  -<  9  <  - 

and 

m  =  mass  of  transfer  vehicle  system 


Figure  2-2.  Thrust  Vector  Location  Relative 
to  the  Orbital  Frame 


system,  the  rate  of  change  in  the  orbital  elements  can  be 
found  utilizing  the  Lagrange  planetary  equations  in  their 
acceleration  component  form  (2:397-407)  : 
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Where 


a^  =  radial  component  of  acceleration 
a^  =  tangential  component  of  acceleration 
a^  =  normal  component  of  acceleration 

Using  the  relation  for  the  distance  from  the  primary 


body  (2:20,24) 


r  = 


1  +  e  cos  v 
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(2-11) 


the  Lagrange  equations  can  be  expressed  in  terms  of  the 
orbital  elements  and  the  true  anomaly  only.  Thus,  substi¬ 
tuting  this  relation,  equations  (2-5)  thru  (2-10)  become 
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The  independent  variable  can  be  changed  from  time,  t, 
to  the  true  anomaly,  v,  using  the  following  relation  for  the 
angular  momentum  (2:17,28) 
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Solving  for  the  time  rate  of  change  in  the  true  anomaly 
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Multiplying  equations  (2-12)  thru  (2-19)  by  results  in 

the  final  form  of  the  perturbation  equations  which  define 
the  changes  in  the  orbital  elements  with  respect  to  v 
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where  the  components  of  acceleration  have  been  expressed  in 
terms  of  the  control  angles. 

Since  the  thrust  is  constant,  and  the  mass  is  assumed 
to  be  constant  over  one  orbit,  the  vehicle  acceleration  may 
be  moved  outside  the  integrals.  In  addition,  with  the 
assumption  that  the  orbital  elements  can  be  treated  as 
constants  over  one  orbit,  they  too  can  be  moved  outside  the 
integrals,  where  possible.  Thus,  the  equations  for  the 
change  in  the  orbital  elements  become 
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Since  this  study  does  not  address  noncoplanar  trans¬ 
fers,  the  out-of-plane  control  angle,  0,  will  be  considered 
to  be  zero  making  equation  (2-34)  for  Ai  and  equation  (2-35) 
for  AO  also  equal  to  zero.  In  addition,  only  the  changes  in 
the  semimajor  axis  and  eccentricity  are  applicable  to  the 
remaining  derivation,  eliminating  the  need  of  equation 
(2-36)  for  A<*>  and  equation  (2-37)  for  AMo. 

Nondimens ional  Form .  The  dependence  of  equations 
(2-32)  and  (2-33)  on  the  semimajor  axis  and  acceleration 
can  be  removed  by  introducing  the  nondimensional  variables 
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Equations  (2-39)  and  (2-40)  represent  the  primary  governing 
equations  for  the  fast  timescale  problem.  Singularities 
occur  in  them  when  the  eccentricity  is  equal  to  zero  or  one, 
e.g.,  if  e  =  0  equation  (2-40)  has  a  singularity.  If  e  =  1, 
both  equations  (2-39)  and  (2-40)  each  have  a  singularity 
when  v  =  n .  However,  only  the  singularity  at  e  =  0  is  of 
concern  since  e  =  1  represents  a  parabolic  trajectory. 

These  singularities  are  removed  in  the  final  algorithm  of 
the  following  section  by  extrapolation  of  the  data  obtained 
from  these  equations . 

Planar  Thrust  Angle  Opt imal  Control  Laws 

Unconstrained  Problem .  Before  discussing  an 
optimal  control  law  for  the  constrained  problem,  we  should 
first  consider  the  unconstrained  problem  to  gain  insight 
into  the  maximum  changes  in  the  semimajor  axis  and  eccentri¬ 
city  that  are  possible  in  the  coplanar  problem.  While  the 
unconstrained  optimal  control  laws  will  provide  useful 
comparison  data,  they  are  obviously  impractical  since  they 
optimize  one  variable  while  allowing  unrestrained  changes  in 
the  others  . 

The  optimal  control  laws  for  the  planar  thrust  control 
angle,  a,  which  maximize  the  change  in  the  semimajor  axis 
and  the  eccentricity  over  one  orbit  in  the  unconstrained 
case  can  be  obtained  using  the  performance  indices  given  by 
(3:4.7-4  8) 
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The  magnitude  of  these  performance  indices  can  be 
maximized  by  first  taking  their  variation  resulting  in 
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where  the  variations  in  the  derivatives  of  the  semimajor 
axis  and  eccentricity  derived  from  equations  (2-20)  and 
(2-21)  are 
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A  stationary  point  (preferably  a  minimum  or  maximum)  is 
obtained  when  the  variance  of  the  performance  index  equals 
zero  (<5J  =  0)  for  all  possible  values  of  the  optimizing 
variable  (Sa) .  This  is  only  true  if  the  integrands  of  equa¬ 
tions  (2-43)  and  (2-44)  are  zero  (3:49).  Thus, 
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Solving  these  two  equations  for  a(y)  yields  the  following 
two  optimal  control  laws  for  the  planar  thrust  angle  in  the 
unconstrained  problem 
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Equation  (2-49)  represents  the  optimal  control  law  for 
maximizing  the  change  in  the  semimajor  axis  and  equation 
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(2-50)  represents  the  optimal  control  law  for  maximizing  the 
change  in  eccentricity.  These  two  control  laws  are  inde¬ 
pendent  of  the  vehicle  acceleration,  the  orbit  semimajor 
axis,  and  the  out  of  plane  thrust  control  angle  (0)  ,  making 
them  valid  for  any  coplanar  or  noncoplanar  transfer.  Direct 
substitution  into  equations  (2-39)  and  (2-40)  will  result  in 
the  maximum  possible  change  in  the  magnitude  of  the  nondi- 
mensional  semimajor  axis  or  eccentricity  over  one  orbit  of  a 
coplanar  transfer. 

Constrained  Problem .  The  fast  timescale  problem 
is  constrained  by  forcing  the  distance  between  the  primary 
focus  (Earth’s  center)  and  the  transfer  orbit  perigee  or 
apogee  position  to  remain  constant.  This  distance  is 
defined  as  (2:25) 
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the  sign  for  the  radius  of  perigee,  rp . 

The  time  rate  of  change  in  r  is  given 
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As  was  done  previously  for  the  orbital  element  perturbation 
equations,  the  independent  variable  can  be  changed  from  t  to 
v  by  multiplying  by 
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Integrating  this  expression  over  one  orbit  results  in 
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Applying  the  assumption  that  the  orbital  elements  and 
acceleration  are  constant,  equation  (2-54)  simplifies  to 

Ar  =  Aa  (1  ±  e)  ±  a  Ae  (2-55) 

<x.  p 


This  expression  can  be  nondimensional ized  in  the  same  manner 
as  done  previously  for  Aa .  Thus,  a  new  nondimensional 
parameter  can  be  defined  by  multiplying  equation  (2-55)  by 
(n2/A}  producing 
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Ar  =  £  Ar  *  J  Aa  (1  ±  e)  ±  Ae  (2-56) 
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By  direct  substitution  of  equation  (2-38)  all  dependence  on 
the  vehicle  acceleration  and  the  orbit  semimajor  axis  can  be 
•  removed  resulting  in  the  final  form  of  the  constraint 

relationship : 
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The  optimal  control  laws  for  a  which  will  maximize  the 
magnitude  of  the  change  in  semimajor  axis  or  eccentricity 
can  be  obtained  by  using  the  performance  indices 
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However,  we  must  now  add  the  constraint  relationship 
described  for  r  ,  making  these  indices  of  the  form  (3:48) 
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where  and  X^  are  the  Lagrange  multipliers.  Substituting 
equation  (2-52)  into  these  two  performance  index  relations 
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Taking  the  variance  of  equations  (2-60)  and  (2-61) ,  noting 

that  Ar  is  constant  (equal  to  zero) ,  and  setting  the 
•  >  p 

integrands  to  zero  (3:49)  yields 
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where  the  variations  in  the  derivatives  of  the  semimajor 
axis  and  eccentricity  are  defined  in  equations  (2-45)  and 
(2-46)  .  Solving  these  equations  for  a(v>)  produces 
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where,  for  the  control  law  a(i>)  which  maximizes  the  change 
in  semimajor  axis,  the  functions  and  Fz  are 
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and,  for  the  control  law  <*(w)^  which  maximizes  the 
eccentricity , 
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By  introducing  a  new  nondimensional  parameter 

X*  =  X  a 

•  • 

equations  (2-67)  and  (2-68)  can  be  rewritten  as 
F1#  =  (l  -  e2)  +  X*  £  2e  (1  ±  e)  ±  (l  -  e2)  J 
F  =  (1  ±  X*)  (l  -  e2) 
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Thus,  with  the  constrained  radii,  two  control  laws  have 
been  found,  one  maximizing  the  change  in  semimajor  axis 
[equations  (2-64)  thru  (2-66)]  and  the  other  maximizing  the 
change  in  eccentricity  [equations  (2-64) ,  (2-69) ,  and 

(2-70)],  that  are  independent  of  the  vehicle  acceleration, 
orbit  semimajor  axis,  and  out-of-plane  thrust  control  angle. 
The  final  solution  to  the  fast  timescale  problem  can  be 
found  by  substituting  these  control  laws  into  equations 
(2-39)  and  (2-40)  to  find  values  of  X  and  X  which  will 

a  e 

*  * 

drive  equation  (2-57)  for  Ar  or  Ar  to  zero.  The 

a  p 

relationship  between  Ar  and  Ar  given  by  (2-56)  shows 

a«p  a,p 

Ar  and  Ar  will  also  be  zero  for  these  values  of  the 

a  p 

Lagrange  multipliers,  thus  satisfying  the  given  constraint 
relationship . 

With  the  constrained  radii  control  laws  derived,  the 
definition  of  even  and  odd  functions  (5:475)  can  be  applied 
to  evaluate  equations  (2-36)  and  (2-37)  for  the  change  in 
the  argument  of  periapsis,  Aco,  and  mean  anomaly  at  epoch, 
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AMo.  The  constrained  control  laws  are  odd  functions  (refer 

to  Figure  2-9) .  With  9  set  to  zero,  insertion  of  any  of  the 

control  laws  into  equation  (2-36)  makes  the  remaining  two 

terms  odd.  Since  the  integral  of  an  odd  function  over  its 

total  period  is  zero,  both  terms  are  zero.  Thus,  over  one 

orbit,  Aw  =  0.  Insertion  into  equation  (2-37)  results  in 

3nt 

the  first  two  terms  being  odd,  leaving  AMo  equal  to  — = —  Aa . 
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Resulting  Algorithms  and  Solutions 

Two  computer  programs  were  written  that  numerically 

*  * 

solved  equation  (2-39)  for  Aa  and  equation  (2-40)  for  Ae  , 

* 

determining  Ar  as  a  function  of  the  Lagrange  multipliers 

a.p 

(for  values  of  eccentricity  ranging  between  zero  and  one) . 

The  first  program  utilized  the  constrained  control  law,  a^, 

* 

given  by  equations  (2-64)  thru  (2-66)  to  maximize  Aa  (Aa  ) . 
The  second  program  utilized  the  const rai  -cntrol  law,  , 

given  by  equations  (2-64) ,  (2-69) ,  and  (2-70)  to  maximize  Ae 

* 

(Ae  ) .  Both  programs  incorporate  the  composite  Simpson's 
rule  (13:156)  to  solve  the  integrals  found  in  the  equations 
for  Aa  and  Ae  .  Thus,  ignoring  the  truncation  error,  the 
equation  used  to  solve  each  integral  is  given  by 

Integral  =  ^Ayj  jl  (u=0)  +  4I(^=Av)  +  2I(u=2Av) 

+  4I(y=3Av)  +  2I(v=4Ai>)  +  ... 

+  2 1  (.v=2  <  tt-Av  >  )  -t-  4 1  (v=2n-Ak>) 

+  I(2tt)\  (2-71) 
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I  =  Integrand 
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The  value  of  n  used  in  the  above  relation  was  set  to  360  in 
each  of  the  runs  of  these  programs  to  allow  the  planar 
thrust  angles  (control  laws)  to  be  calculated  at  one  degree 
intervals  around  the  orbit  .  Higher  values  of  n  did  not 
improve  the  accuracy  within  the  limits  of  the  machine  used 
<5=  1  .E-16)  . 


These  two  programs  were  used  to  obtain  a  rough  estimate 

♦  <  *  ) 

of  the  range  of  the  Lagrange  multipliers,  X^  and  X (V  ) 

*  *  * 

where  the  the  functions  Ar  and  Ar  cross  the  X.  axis  (equal 

Cl  p  l 

to  zero) .  Sample  data  (plotted  in  Figures  2-3  thru  2-6) 
shows  that  the  desired  values  of  these  parameters  lie 
between  ±  1.  This  range  provided  the  initial  guesses  to  the 
secant  method  routine  of  Program  DELAMAX2  (Appendix  B-l)  and 
Program  DELEMAX2  (Appendix  B-2) ,  expansions  of  the  original 
two  programs.  Given  two  initial  guesses  in  the  vicinity  of 
the  solution  the  next  guess  of  X*  >  is  found  from  (13:265) 


<*> 


j  ♦  * 


=  X 


<  *  > 


Ar 


a.P 


X<*>  -  x( 


j-i 


Ar 


a,p 


-  Ar 


a,p 


(2-72) 


j-* 


where  the  subscript  j  represents  the  current  value  of  the 
parameter.  Thus,  with  this  routine,  these  two  programs 
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Figure  2-3. 


Nondimensional  Change  in  Radius  of  Apogee  as 
a  Function  of  the  Lagrange  Multiplier  (max 
Aa  Control  Law)  at  Various  Eccentricity 
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Figure  2-4  , 


Nondimensional  Change  in  Radius  of  Perigee  as 
a  Function  of  the  Lagrange  Multiplier  (max 
Aa  Control  Law)  at  Various  Eccentricity 
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Figure  2-5 . 


Nondimensional  Change  in  Radius  of  Apogee  as 
a  Function  of  the  Lagrange  Multiplier  (max 
Ae  Control  Law)  at  Various  Eccentricity 
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Figure  2-6, 


Nondimensional  Change  in  Radius  of  Perigee  as 
a  Function  of  the  Lagrange  Multiplier  (max 
Ae  Control  Law)  at  Various  Eccentricity 


solved  the  fast  timescale  problem,  providing  X.^,  ,  Aa  , 

* 

and  Ae  as  functions  of  e,  while  holding  or  Arp  equal  to 

zero.  Resulting  plots  of  these  parameters  are  provided  in 
Figures  2-7  thru  2-9. 

The  uniqueness  of  the  particular  constraints  imposed  is 

0  0  0  0 

seen  in  the  plots  of  Aa^,  ±Aefl ,  ±Aa^  ,  and  Ae^  displayed  in 
Figures  2-8  and  2-9.  Notice  that  imposing  each  of  the  two 
constraints  results  in  identical  changes  in  the  magnitudes 

Ik 

of  Aa  and  Ae  ,  whether  these  changes  were  obtained  from  the 

* 

control  law  a  which  maximizes  Aa  or  the  control  law  a 

a  e 

which  maximizes  Ae  .  However,  while  the  magnitudes  of  these 
resulting  changes  are  the  same  for  each  of  the  two  control 


e 

Figure  2-7,  Lagrange  Multiplier  as  a  Function  of 
Eccentricity  for  Constant  Apogee 
and  Perigee 
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Figure  2-8.  Nondimensional  Change  in  Semimajor  Axis 
as  a  Function  of  Eccentricity  for 
Constant  Apogee  and  Perigee 
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Figure  2-9.  Change  in  Eccentricity  Parameter  as  a 
Function  of  Eccenticity  for  Constant 
Apogee  and  Perigee 


laws,  when  constrained  by  Ar  =  0,  a  produces  changes  of 

Cl  a 

opposite  sign  (negative)  to  those  found  using  . 

Two  similar  programs  incorporating  the  composite 
Simpson's  rule  algorithm,  DELAMAXUC  (Appendix  B-3)  and 
DELEMAXUC  (Appendix  B-4) ,  provide  the  solutions  to  equations 
(2-39)  and  (2-40)  for  the  unconstrained  problem  utilizing 
the  control  laws  defined  in  equations  (2-49)  and  (2-50) . 

The  first  program  (utilizing  a  e)  provides  the  maximum 

* 

possible  changes  in  the  nondimensional  semimajor  axis,  Aa^, 
and  the  associated  residual  changes  in  the  eccentricity 

4|t 

parameter,  Ae  ,  as  functions  of  eccentricity.  The  second 

Cl 

(utilizing  a  )  provides  the  maximum  possible  values  of  the 

euc  4 

m 

changes  in  the  eccentricity  parameter,  Ae^,  and  the 
associated  changes  in  the  nondimensional  semimajor  axis, 

Aa^ ,  as  functions  of  eccentricity. 

The  data  obtained  from  these  two  programs  is  plotted  in 
Figures  2-8  and  2-9  for  comparison  with  the  constrained 
results.  The  plots  show  the  maximum  obtainable  changes  in 
the  nondimensional  semimajor  axis  or  eccentricity  parameter 
are  reduced  as  eccentricity  is  increased.  In  addition,  it 
can  be  seen  that  utilizing  a  to  increase  the  semimajor 

CLUC 

* 

axis  of  eccentric  orbits  results  in  negative  values  in  Ae^ 

driving  the  eccentricity  to  zero.  Once  the  orbit  is 

* 

circularized,  Ae  becomes  zero  and  the  maximum  change  in  the 

a. 

* 

nondimensional  semimajor  axis  (Aa  =  4rr)  is  obtained.  This 

a 

maximum  value  is  identical  to  the  result  obtained  by  Alfano 


for  spiral  transfers  between  coplanar  circular  orbits  where 
both  a  and  6  are  constant  and  equal  to  zero.  Thus,  this 
control  law  is  ideal  to  increase  the  semimajor  axis  while 
circularizing  an  eccentric  orbit.  There  is  never  a  threat 
of  impact  with  the  primary  body,  since  the  changes  in  Ar 

Cl 

and  Ar^  are  both  positive.  On  the  other  hand,  utilizing 
this  control  law  to  decrease  the  semimajor  axis  (thrusting 
in  the  opposite  direction)  results  in  an  uncontrollable 
increase  in  the  eccentricity  and  decrease  in  Ar  and  Ar  , 

a.  p 

making  it  impractical . 

Utilizing  a  to  increase  the  eccentricity  of  a  cir- 
cular  or  eccentric  orbit  produces  positive  changes  in  the 
semimajor  axis.  However  Ar^  is  negative,  leading  to  a 
possible  impact  with  the  primary  body.  Thrusting  in  the 
opposite  direction  decreases  both  the  eccentricity  and 
semimajor  axis  of  an  eccentric  orbit.  However,  Ar  is 

CL 

negative  while  Arp  is  positive,  making  an  application  to 
circularizing  eccentric  orbits  possible  if  the  uncontrol¬ 
lable  changes  in  the  semimajor  axis  were  acceptable. 

Unfortunately,  the  constraint  relationship  imposed  in 
this  study  provides  results  significantly  less  than  the 
optimal.  Figure  2-8  shows  when  e  =  0,  the  magnitude  of 
change  in  the  semimajor  axis  is  slightly  higher  than  one 

half  (=  2.2tt)  of  that  obtained  in  the  unconstrained  case. 

* 

As  e  — *  1 ,  the  magnitude  of  Aa  constrained  by  Ar^  =  0 
approaches  the  unconstrained  value  while  the  magnitude  of 


m 

Aa  constrained  by  Ar  =0  decreases  to  zero.  Figure  2-9 

a 

0 

shows  similar  results  for  Ae .  When  e  =  0,  the  magnitude  of 
change  in  the  eccentricity  is  exactly  the  same  as  the  change 
in  semimajor  axis  and  approximately  two  thirds  of  the 
unconstrained  value.  As  e  — ►  1,  the  magnitude  of  Ae* 
constrained  by  Ar^  =  0  approaches  the  unconstrained  value 
while  the  magnitude  of  Aa  constrained  by  Ar^  =  0  diverges 
toward  zero. 

The  data  from  Figures  2-8  and  2-9  can  be  used  in 
conjunction  with  the  vehicle  accelerations  derived  in 
Appendix  A  to  validate  that  the  changes  in  eccentricity  and 

semimajor  axis  over  one  orbit  are  small.  Utilizing  the 

0  0 

maximum  values  of  Aa  and  Ae  for  the  unconstrained  control 
laws  (4n  and  S*  3 .  In ,  respectively)  and  the  maximum  value  of 
acceleration  seen  by  the  arcjet  TVS  at  low  earth  orbit  (LEO) 
and  geosynchronous  orbit  (GEO) ,  equation  (2-38)  for  the 
dimensional  values  of  Aa  and  Ae  can  now  be  solved: 


Aa  = 


A  A  * 
—  Aa 
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A  a  * 


A  a' 


3.986  x  1014  m3/sZ 


(4rc) 


*  (3.153  x  10"14  m-2)  A 


Ae  = 


A  .  * 
Ae 
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n  a 


A  a  * 
-  Ae  = 


A  a 


3.986  x  1014  m3/s2 


(3.  ITT) 


=  (2.443  x  10"14  m~2} 


A  a2 


At  LEO  (a  =  6.68  x  lO^m.  A  =  ,00446  m/s2)  these  equations 

max 


yield  a  change  in  semimajor  axis  of  =  .6%  of  its  original 


value  and  a  total  change  in  eccentricity  is  =  .005.  At  GEO 
(a  =  4.22  x  10?m,  A  =  .00257  m/s2),  the  change  in 

TflOX 

semimajor  axis  has  increased  to  =  14.595  of  its  original 
value  and  the  total  change  in  eccentricity  is  up  to  .112. 
Thus,  as  expected,  the  original  assumption  is  valid  near 
LEO,  but  becomes  much  less  accurate  as  the  semimajor  axis 
increases  to  GEO. 

The  programs  of  Appendices  B-l  and  B-2  were  also  used 
to  obtain  the  changes  in  the  constrained  planar  thrust  angle 
(control  laws) ,  or  ,  as  the  spacecraft  proceeds  around  an 
orbit  of  given  eccentricity.  Data  for  orbits  with  eccentri¬ 
cities  of  .01,  .4,  and  .8  is  plotted  in  Figure  2-10. 


Figure  2-10,  Constrained  Thrust  Vector  Angle  as  a 
Function  of  the  True  Anomaly  for 
Various  Eccentricity 


For  the  constraint  relationship  Arp  =  0,  the  plots  of 
a  and  a  only  differ  slightly  at  e  =  .01  and  become 

<X  # 

indistinguishable  as  the  eccentricity  becomes  larger.  When 
apogee  is  constrained  (Ar  =  0) ,  the  plots  of  a  and  a 
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differ  by  a  shift  of  ±180  degrees  (again,  slightly  more  for 
lower  values  of  e) .  This  "phase  shift"  explains  the 
difference  in  signs  of  Aa  and  Ae  obtained  by  these  two 
control  laws  earlier  (see  Figures  2-8  and  2-9)  when  apogee 
was  constrained.  By  simply  shifting  the  curves  for  one  of 
the  control  laws  by  180  to  more  closely  agree  with  the 

curves  of  the  other  we  can  obtain  the  change  in  sign  of  Aa 

* 

and  Ae  needed  to  make  the  results  match.  Thus,  even  though 
the  two  control  laws  differ  slightly  at  lower  eccentri¬ 
cities,  they  do  in  fact  produce  identical  results  for  the 
changes  in  these  two  parameters  (magnitude  and  sign) . 

A  final  point  regarding  the  plots  of  Aa  and  Ae  in 
Figures  2-8  and  2-9  involves  their  significance  when  related 
to  a  desired  transfer.  For  example,  to  increase  the  radius 
of  a  circular  orbit  (increase  a)  using  the  specified 
constraints,  perigee  is  first  held  constant  while  pushing 
apogee  out  to  the  larger  radius.  This  would  result  in  an 
increase  in  eccentricity  since  we  are  basically  "stretching" 
the  original  orbit  from  a  circular  to  an  elliptical  shape. 

To  complete  the  transfer,  apogee  is  held  constant  while 
raising  perigee  to  the  final  radius,  decreasing  the 
eccentricity  until  the  orbit  was  circularized  (e  =  0) .  On 
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the  other  hand,  to  decrease  the  radius  of  a  circular  orbit 
(decrease  a),  this  procedure  is  reversed.  First,  apogee  is 
held  constant  and  perigee  is  lowered  to  the  desired  radius, 
increasing  the  eccentricity.  Then,  perigee  is  held  constant 
while  apogee  is  pushed  down  to  the  final  radius,  decreasing 
the  eccentricity  to  zero.  Therefore,  it  is  expected  that 
the  control  laws  constrained  by  maintaining  perigee  constant 
will  produce  changes  in  the  parameters  Aa*  and  Ae*  which  are 

of  the  same  sign  (±) ,  while  the  control  laws  constraining 

*  _  * 

apogee  produce  changes  of  the  opposite  sign  (±Aa  and  +Ae  )  . 
Figures  2-8  and  2-9  show  this  is  exactly  what  has  been 

found.  For  both  control  laws,  constraining  perigee  results 

*  * 

xn  positive  changes  in  Aa  and  Ae  .  If  is  put  "in  phase" 

with  a* ,  constraining  apogee  produces  positive  changes  in 

*  * 

Aa  and  negative  changes  in  Ae  for  both  control  laws. 

o 

Thrusting  in  the  opposite  direction  (shifting  or  ±180  ) 
changes  the  sign  of  all  of  the  parameters  making  both  Aa 

4k 

and  Ae  negative  when  constraining  perigee,  and  Aa  negative 
and  Ae  positive  when  constraining  apogee. 

To  better  illustrate  the  final  four  control  laws  (or 

V 

o 

and  or  ±  180  for  each  of  the  two  constraints) ,  Figures 
2-11  (a)  and  2-11  (b)  provide  sketches  of  the  required 

thrust  vector  control  for  a  single  orbit  utilizing  the  data 
(e  =  .40)  provided  in  Figure  2-10.  When  perigee  is 
constrained  [Figure  2-11  (a)],  the  direction  of  the  thrilst 
vector  changes  very  rapidly  near  apogee.  When  apogee  is 
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Constant  Perigee 


(b)  . 


Constant  Apogee 


Figure  2-11. 
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A> 


constrained  [Figure  2-11  (b) ] ,  the  rapid  change  occurs  near 
perigee.  These  characteristic  can  also  be  seen  in  the 
slopes  of  the  plots  of  ci^  in  Figure  2-10. 

To  complete  the  short  timescale  problem,  a  final 
program,  Program  INTERPO  (Appendix  B-5) ,  was  written  to 
provide  a  solution  for  che  A*  \  Aa*  and  Ae*  for  any  given 
value  of  eccentricity  less  than  one.  This  includes 
solutions  at  e  =  0  that  were  unobtainable  in  the  previous 
programs  due  to  the  singularity  in  equation  (2-40) .  A 
Newton  formula  (13:92-98)  was  incorporated  to  produce  an 
interpolating  or  extrapolating  nth  order  polynomial  to 
calculate  the  values  of  the  parameters  between  the  data 
points  obtained  from  the  earlier  programs  which  solved 
equations  (2-39)  and  (2-40) .  This  program  is  utilized  as  a 

subroutine  in  the  slow  timescale  solution  algorithm  to 

*  0 

provide  an  efficient  method  of  calculating  Aa  and  Ae  as 
functions  of  eccentricity. 


III.  The  Slow  Timescale  Problem 


Problem  Statement 

The  slow  timescale  problem  provides  for  the  complete 
transfer  of  a  spacecraft  from  an  initial  tc  final  orbit 
through  many  revolutions  about  the  primary  body.  Changes  in 
the  orbital  elements  found  in  the  fast  timescale  problem  of 
chapter  2  are  included,  as  well  as  the  changes  in  the 
vehicle  mass  due  to  propellant  expulsion. 


Derivat i on 

Perturbat i on  Equations  ■  The  period  of  an  elliptical 
orbit  is  found  from  <2:33) 

TT£P  =  —  (3-1) 


Given  the  low  thrust  assumption,  the  changes  in  the  semi¬ 
major  axis  and  eccentricity  with  respect  to  time  can  be 
approximated  from  the  single  orbit  case  of  the  fast  time- 
scale  solution  by 
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where  the  nondimensional  parameters,  Aa  and  Ae  , 


are  given 


in  equations  (2-39)  and  (2-40) .  The  "sign"  of  each  equation 


is  dependant  on  the  change  in  the  unconstrained  radii 
desired  (increase  or  decrease)  and  the  control  law  chosen 
(a  or  a  )  .  To  remove  the  dependence  of  the  above  two 

a  • 

equations  on  the  vehicle  acceleration,  the  independent 
variable  can  be  changed  from  time  to  velocity: 
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Based  on  the  initial  semimajor  axis,  a  ,  and  initial  mean 
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anomaly , 

"o' 

we 

can 

define  three 

nondimensional  variables: 

where 
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Substituting  these  parameters  into  equations  (3-4)  and  (3-5) 
produces  the  final  form  of  the  two  governing  equations: 
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(3-9) 


In  addition,  substitution  into  equation  (2-51)  for  the 
radius  of  apogee  and  perigee  results  in  a  new  nondimens i onal 
form  of  the  radii  equation: 


r 

r  =  =  a  (  1  ±  e)  (3-10) 

a,p  a 

o 

Mass  FI ow  Rate/Acceleration  Relation.  The  vehicle  mass 
and  acceleration  are  related  through  Newton's  second  law  of 
mot  ion  (2 : 3-4)  : 


where 


T  =  m  (t)  A  (t) 

T  =  vehicle  thrust  (constant) 
m(t)  =  vehicle  mass  at  time  t 
A(t)  =  vehicle  acceleration  at  time  t 


(3-11) 


At  the  beginning  of  the  transfer  (t  =  t  ,  .  .  =  0)  , 

equation  (3-11)  can  be  written  in  terms  of  the  initial 
vehicle  mass  and  acceleration: 

T  =  m  A  (3-12) 

O  O 

Additionally,  the  vehicle  mass  at  time  t  can  be  expressed  in 
terms  of  the  initial  mass  and  propellant  mass  flow  rate, 
m,  using  the  relation 

m (t)  =  m  -  m  t  (3-13) 

O 


3-3 


Assuming  constant  thrust  implies  the  mass  flow  rate  is  also 
constant  . 


Combining  equations  (3-11)  thru  (3-13)  the  vehicle 
acceleration  can  now  be  modeled  in  the  slow  timescale 
problem  as 

,  A 

A  (t)  =  =  - 2 -  (3-14) 

dt  (i  -  lit) 


where  M  is  the  specific  mass  flow  rate,  m/m  . 

O 

Total  Transfer  Velocity  Change/Time  Relation .  The 
total  accumulated  velocity  change  is  found  by  integrating 
equation  (3-14)  between  0  and  t  resulting  in 

A 

v(t  )  -  v  =  Av (t  )  - - -  In  (1  -  Mt,)  (3-15) 

f  °  f  li  f 


Solving  this  expression  fcr  the  total  transfer  time  provides 
the  final  equation 
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Resulting  Algorithms 

Slow  T imescal e  Constraint  AppI icat i on .  One  method  of 
completing  the  transfer  of  the  slow  timescale  problem 
involves  applying  the  results  of  the  fast  timescale  problem 
for  each  constraint  over  an  extended  period  to  relocate  the 
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perigee  and  apogee  positions  separately.  Program  TRANSMUL, 
Appendix  C-l,  was  written  to  perform  the  complete  transfer 
by  constraining  one  radii  (apogee  or  perigee)  over  many 
revolutions  until  the  other  radii  matches  the  final  orbit. 
For  example,  if  final  perigee  is  higher  than  that  of  the 
original  orbit,  it  is  first  held  constant  while  apogee  is 
increased  or  decreased,  as  required,  to  its  final  value. 

The  transfer  is  completed  by  constraining  apogee  as  perigee 
is  raised  to  its  final  position.  On  the  other  hand,  if 
final  perigee  is  lower,  apogee  is  first  constrained  while 
perigee  is  lowered,  then  perigee  is  constrained  while  apogee 
is  raised  or  lowered  to  complete  the  transfer.  This 
procedure  applies  to  any  transfer,  including  circular-to- 
circular  transfers.  Thus,  given  the  initial  and  final  orbit 
semimajor  axes  and  eccentricities,  primary  body  gravita¬ 
tional  parameter  and  radius,  and  vehicle  parameters 
(specific  impulse,  initial  mass,  and  propellant  mass  flow 
rate) ,  this  program  calculates  the  total  accumulated 
velocity  change  and  time  for  the  constrained  radii  transfer. 

The  program  incorporates  an  ordinary  differential 
equations  integrator  equipped  with  a  fourth  order  predictor- 
corrector  algorithm,  to  numerically  solve  equations  (3-8) 
and  (3-9) .  In  addition,  as  stated  in  chapter  two,  the  final 
algorithm  of  the  fast  timescale  problem  (Appendix  B-5)  has 
been  incorporated  to  provide  nondimens i onal  changes  in  the 
orbital  elements  needed  to  solve  equations  (3-8)  and  (3-9) . 
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Given  the  current  value  of  the  orbit  eccentricity  and  the 

constraint  condition  (either  constant  r  or  r  ) ,  values  of 

a  p 

4»  4 1 

Aa  and  Ae  are  generated  from  data  obtained  using  the 

constrained  control  law  that  maximizes  Ae  (a  ) .  The  other 

constrained  control  law  that  maximizes  Aa  (ot  )  could  have 

a 

been  used  since  it  provides  identical  results. 

As  a  comparison  in  circular-to-circular  transfers, 
this  program  also  provides  the  Av  requirements  of  spiral 
and  Hohmann  transfers.  For  planar  spiral  transfers,  the 
velocity  change  is  given  simply  by  the  difference  between 
the  velocities  associated  with  each  orbit  (10:463): 


Av  ,  =  I  v  -  v  ! 

spiral  1  csi  cs2 


(3-17) 


where  (2:165-166) 


The  subscripts  (1)  and  (2)  in  the  above  relation  refers  to 
the  initial  and  final  orbit,  respectively.  For  Hohmann 
transfers,  the  total  velocity  change  is  obtained  from 
(2:163-166) 


where 


Av  ,  =  Av  4-  Av 

Hohmann  1  2 


Av.  =  v.  -  v 

i  i  csi 


(3-18) 


v  =  2m 


"  N 

n  1  1 

2u  -  -  - 

®  r  r  +  r 

L  l  «■  1  2  J 


and  r  is  the  orbit  radius 


3-6 


Fast  T imescale  Constraint  Application .  Figures  2-8  and 


2-9  show  when  the  eccentricity  is  zero,  the  changes  in  the 
nondimens ional  semimajor  axes  resulting  from  each  constraint 
are  equal  in  magnitude  and  sign  while  the  changes  in  the 
eccentricity  parameters  are  equal  in  magnitude,  but  opposite 
in  sign.  This  means  that  if  we  apply  one  constraint  for  one 
orbit,  then  the  other  constraint  during  the  next  orbit,  over 
two  revolutions,  the  eccentricity  can  be  held  nearly 
constant  (=  0)  as  the  semimajor  axis  is  increased  or 
decreased.  Program  TRANSALT,  Appendix  C-2 ,  utilizes  this 
idea  to  complete  the  slow  timescale  transfer  by  applying  the 
fast  timescale  results  to  two  revolutions  of  the  transfer  at 
a  time,  instead  of  one  as  was  done  in  the  previous  program. 
Perigee  is  constrained  during  one  of  the  two  revolutions  and 
apogee  during  the  other.  All  routines  of  the  previous 
program  have  been  incorporated,  as  well  as  the  following 
modified  versions  of  equations  (4-8)  and  (4-9)  to  evaluate 
the  average  change  in  the  orbital  elements  occurring  during 
the  two  orbits: 


da  __ 

3/2 

(a) 

r  ±  Aa*  *■ 
p 

Aa* 

a 

dv 

2  n 

2 

de 

1/2 

(I) 

f  4  K  * 

* 

Ae 

a 

dv 

2tt 

l  2 

(3-19) 


(3-20) 


where  the  subscripts  (a)  and  (p)  refer  to  constrained 
perigee  or  apogee,  respectively.  The  "signs"  in  these 


relations  assume  the  control  law  a  is  utilized.  The  upper 
signs  are  used  for  transfers  to  higher  orbits,  while  the 
lower  signs  are  used  for  transfers  to  lower  orbits.  This 
combination  of  constraints  keeps  each  of  the  revolutions 
nearly  circular  as  the  orbit  radius  is  raised  or  lowered, 
making  this  routine  valid  for  circular-t o-circular  transfers 
only  . 

As  an  example,  consider  a  transfer  to  a  higher  orbit. 
During  the  first  revolution  of  the  transfer,  perigee  is 
constrained  while  apogee  is  raised.  At  the  beginning  of  the 
second  revolution,  the  constraint  is  switched  making  apogee 
constrained  while  perigee  is  raised.  By  the  end  of  the 
second  revolution,  the  orbit  has  been  recircularized  at  a 
higher  radius.  This  procedure  is  repeated  until  the  radius 
matches  that  of  the  fina^  desired  orbit. 
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Circular-t o-Circular  Transfers 


Computer  runs  were  performed  for  ratios  of  the  final  to 
initial  orbit  semimajor  axes  between  0  and  100.  Utilizing 
the  data  obtained,  a  comparison  among  the  nondimensional  Av 
requirements  for  each  of  the  constrained  radii,  spiral,  and 
Hohmann  transfers  is  shown  in  Figure  4-1  below.  Sample  data 
for  the  LEO  to  GEO  transfer  performed  using  the  slow  and 
fast  timescale  algorithms  is  provided  in  Appendices  D-l  and 
D-2 ,  respectively.  As  expected  after  reviewing  the  results 
of  the  fast  timescale  problem,  the  Av  requirements  for  the 


Qf/ao 


Figure  4-1  Nondimensional  Total  Accumulated  Velocity 
Change  for  Ratios  of  Final  to  Initial 
Semimajor  Axis 


constrained  radii  transfers  are  higher  than  those  of  the 
spiral  transfer. 

The  plots  of  Figure  4-1  show  the  Av  requirements 
resulting  from  applying  each  constraint  separately  during 
the  slow  timescale  problem  are  slightly  lower  than  those 
obtained  by  applying  a  combination  of  the  constraints  in  the 
fast  timescale  problem.  The  reason  for  this  difference  can 
be  seen  by  referring  to  the  plots  of  Aa*  in  Figure  2-8. 

While  the  fast  timescale  application  keeps  the  eccentricity 
near  zero,  the  slow  timescale  application  allows  the  eccen¬ 
tricity  tc  increase.  Figure  2-8  shows  higher  eccentricity 
results  in  larger  values  of  Aa*  when  perigee  is  constrained 
and  lower  values  of  Aa  when  apogee  is  constrained.  How¬ 
ever,  the  data  obtained  using  the  slow  timescale  application 
show  a  apogee  is  constrained  less  than  30%  of  the  total 
transfer  time.  This  means  that  during  nearly  70%  of  the 
total  transfer,  the  change  in  semimajor  axis  per  revolution 
is  greater  than  that  obtained  from  the  transfer  using  the 
fast  timescale  application  where  the  eccentricity  is 
continually  kept  near  zero.  Thus,  the  application  of  the 
constraints  separately  in  the  slow  timescale  problem  result 
in  a  quicker  and  less  expensive  transfer. 

Figure  4-2  shows  the  resulting  trajectory  of  a  transfer 
between  LEO  and  GEO  using  the  slow  timescale  application  of 
the  constrained  radii  control  law.  Unfortunately,  this 


Circular-to-Eccentric  Trans f ers 


The  most  feasible  application  of  the  constrained  radii 
transfer  is  to  transfers  between  circular  and  eccentric 
orbits.  One  such  transfer  is  that  needed  by  communications 
satellites  designed  to  provide  coverage  of  the  northern 
hemisphere  (9:54.3.1-54.4.5).  These  satellites  are  placed 
into  Molniya  orbits,  named  after  the  Soviet  communications 
satellites  that  first  used  them  in  the  60' s.  Molniya  orbits 
are  highly  elliptical  (e  =  .73,  depending  on  application) 
with  a  period  of  12  hr  and  an  inclination  of  63.4  degrees. 

The  inclination  of  Molniya  orbits  is  specifically 
chosen  to  maintain  the  argument  of  periapsis,  co,  at  270 
degrees,  keeping  apogee  directly  above  the  northern 
hemisphere.  The  major  perturbation  influencing  satellites 
in  earth  orbit  is  due  to  the  earth’s  oblateness,  the  term 
in  the  geopotential  (14:46,86-91).  By  substituting  the 
secular  terms  from  the  Jz  disturbing  function  into  the 
disturbing  function  form  of  the  Lagrange  planetary  equation 
for  w 


to  = 

it  can  be  seen 
the  right  hand 
this  critical 


3  n  J. 


2  a  (1 


If_  r 

2  2  I 
-  e  )  v 


5.2.  0 

sm  i-2 


(4-1) 


that  by  choosing  an  inclination  of  63,4  , 
term  becomes  zero,  driving  u>  to  zero.  Thus, 
inclination  eliminates  the  influence  of  the 


earth's  oblateness  on  the  spacecraft. 


# 


A  simulated  transfer  to  a  Molniya  orbit  was  performed 
using  the  5000  kg  arcjet  transfer  vehicle  (discussed  in 
Appendix  A) .  The  transfer  assumed  an  initial  807  km  circular 
orbit  (a  =  7185  km)  with  inclination  of  63.4°.  Thrusting 
begins  at  u  =  270°  using  the  constrained  perigee  control  law 
(w  =  0,  initially)  and  continues  until  the  semimajor  axis 
and  eccentricity  are  increased  to  the  final  Molniya  orbit 
(a  =  26,610  km  and  e  =  .73) .  Data  obtained  from  the  comput¬ 
er  run  is  included  in  Appendix  D-3 .  In  addition,  Figure  4-3 
below  provides  a  plot  of  the  spacecraft  trajectory  obtained 
to  illustrate  the  transfer  performed. 

The  data  obtained  indicates  the  transfer  to  the  Molniya 
orbit  required  a  total  Av  of  5.83  km/sec  and  took  approxi- 


Figure  4-3.  Spacecraft  Trajectory  for  Transfer 
to  Molniya  Orbit  using  Constrained 
Radii  Control  Law 


Note  : 


The  revolution  number  is  indicated  where  possible. 
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mately  72.4  days  to  complete.  With  these  values  known, 
using  equation  (A-15)  the  total  propellant  used  is 


m  =  mt  =  (.398  x  10-3  kg/sec) (3.13  x  10*  s)  =  1244  kg 

p 

For  comparison,  consider  executing  this  transfer  with 
chemical  system.  Only  the  first  half  of  the  Hohmann 
transfer  needs  to  be  performed,  thus  from  equation  (3-18) 
the  total  Av  is  found  by 

Av  =  I  v  -  v 

1  1  C3l 


where  r  and  r  are  the  final  orbit  radius  of  apogee  and 
a  p 

perigee,  respectively.  This  equation  yields  a  Av  of  2.35 
km/sec.  The  total  propellant  needed  can  be  found  using 
equation  (A-7) 

1  -  e*P  1 

v  ■P  C''  J 

Thus,  assuming  an  I  of  300  sec,  this  equation  gives  a 
propellant  mass  of  2749  kg,  over  twice  that  needed  by  the 
arcjet  propelled  vehicle. 


Eccentric-to-Eccentric  Trans f er s 


While  the  slow  timescale  constraint  application 
algorithm  was  designed  for  application  to  eccentric-to- 
eccentric  orbit  transfers,  it  provides  no  control  of  <*>  other 
than  keeping  it  constant.  Thus,  it  is  only  useful  for 
eccentric  transfers  which  require  no  change  in  w.  The 
resulting  trajectory  for  such  a  transfer  would  be  similar  to 
the  circular-t o-circular  transfer  shown  in  Figure  4-2. 


Cone lus i ons  and  Recommendations 


V  . 

Two  optimal  control  laws  for  a  continuous,  low  thrust 
spacecraft,  each  resulting  in  identical  changes  in  the 
orbital  elements  have  been  derived.  Constrained  by  constant 
radius  of  perigee  or  apogee,  these  control  laws  can  be 
applied  to  provide  an  optimal  coplanar  circular-t o-eccentric 
or  eccentric-to-eccentric  (Aco  =  0)  transfer,  supplying 
needed  control  of  the  perigee  and  apogee  heights.  Appli¬ 
cation  to  circular-to-circular  transfers  is  not  feasible 
since  it  is  much  more  expensive  (larger  Av)  and  complex 
[a  =  a  (l>)  ]  to  perform  than  the  optimal  spiral  control  law 
(«  =  0)  . 

The  results  of  the  application  of  the  constrained  radii 
control  laws  in  the  slow  timescale  problem  inherently  become 
less  accurate  as  the  semimajor  axis  increases.  This  is  due 
to  the  associated  increase  in  orbital  period  which  slowly 
makes  the  change  in  semimajor  axis  more  and  more  signifi¬ 
cant.  Therefore,  the  assumption  of  the  fast  timescale  which 
states  that  the  changes  in  the  orbital  elements  are  small 
eventually  becomes  invalid.  However,  as  was  the  case  with 
the  control  law  derived  by  Alfano  (1:52),  in  actual  appli¬ 
cation,  a  closed  loop  guidance  scheme  would  be  implemented 
during  the  latter  part  of  the  transfer  using  actual  values 
of  the  orbital  elements.  This  would  also  eliminate  any 

accumulated  during  the  earlier  part  of  the  transfer. 


errors 


A  recommendation  to  further  this  study  is  to  add  the 
change  in  the  argument  of  periapsis,  Aw,  to  the  constraint 
relationship.  This  would  allow  for  control  of  <*>  during 
transfers  between  planar  eccentric  orbits  that  is  not 
included  by  this  analysis.  In  addition,  other  approaches  t 
optimizing  the  eccentric  transfer  problem  should  be  consid¬ 
ered  in  hopes  of  achieving  better  performance.  For  example 
the  development  of  an  optimal  control  law  to  provide 
specified  changes  in  perigee,  apogee,  and  the  argument  of 
periapsis.  However,  successful  completion  of  such 
analyses  would  depend  on  the  complexity  of  the  slow 
timescale  problem. 


Appendix  A:  Perf ormance  Analysis  of  Proposed  Electrical 
Propul s i on  Systems  for  an  Earth  Orbiting  Trans f er  Vehicle 


The  following  analysis  utilizes  information  provided  in 
a  paper  presented  to  the  1985  Joint  Army-Navy-NASA-Air  Force 
Propulsion  Meeting  by  Smith  and  Knowles  discussing  studies 
of  electrical  propulsion  systems  (10:457,463-467).  Supple¬ 
mental  information  on  systems  currently  in  development  was 
provided  by  Rocket  Research  Company  and  two  leading  electric 
propulsion  research  centers:  NASA  Lewis  Research  Center  for 

ion  systems  and  the  Air  Force  Astronautics  Laboratory  for 
arcjet  systems. 

The  purpose  of  this  analysis  is  to  estimate  the  maximum 
acceleration  placed  on  the  transfer  vehicle  system  (TVS)  and 
the  maximum  percent  change  in  mass  of  the  TVS  over  one 
orbit.  The  analysis  considers  a  typical  transfer  vehicle 
mission  of  delivering  a  payload  from  low  earth  orbit  (LEO 
300  km)  to  geosynchronous  orbit  (GEO  ~  35,863  km) ,  then 
returning  to  LEO  in  preparation  for  the  next  mission.  For 
comparison,  two  types  of  electrical  propulsion  systems  are 
considered  for  use  by  the  transfer  vehicle:  an  Ammonia 

(N  H  )  arcjet  system  and  a  Xenon  (Xe)  ion  system. 

2  ■* 

The  Solar  Array  Flight  Experiment  (SAFE)  flown  on 
STS-41D  has  demonstrated  that  a  solar  power  system  can  be 
built  to  provide  a  mass-to-power  ratio  or  approximately  15 
kg/kW  (7) .  Based  on  current  technology  (8;  12) ,  the  arcjet 
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propulsion  system  could  be  equipped  with  one  27  kW  N2H^ 
thruster  having  a  mass  of  approximately  4  kg,  specific 
impulse  (I  ) ,  of  863  sec,  and  an  efficiency  of  37%.  With 
the  given  performance  parameters,  this  thruster  delivers 
3.367  N  of  thrust  CIO :456 .Table  IV) .  It  may  be  possible  for 
the  power  conditioning  unit  (PCU)  associated  with  this 
thruster  to  have  an  efficiency  of  95%  and  a  mass  as  little 
as  15  kg.  Cooling  of  the  PCU  could  be  performed  by  heat 
exchange  with  the  propellant,  eliminating  the  need  for 
radiators.  Thus,  assuming  10%  of  the  available  power  is 
needed  for  supporting  electronics  and  losses  in  the  PCU,  the 
total  power  requirement  for  the  arcjet  transfer  vehicle  is 
30  kW .  This  makes  the  mass  of  the  supporting  power  system 
approximately  450  kg.  Allowing  an  additional  6  kg  for 
connecting  hardware,  the  total  propulsion  system  mass  is  25 
kg  and  the  dry  mass  of  arcjet  transfer  vehicle,  excluding 
fuel  tanks  and  structure,  becomes  475  kg. 

For  the  transfer  vehicle  equipped  with  a  Xe  ion 
propulsion  system,  five  10.87  KW  thrusters  will  be 
considered,  each  with  an  I  of  4267  sec  and  an  efficiency 

SP 

of  75%,  providing  389  mN  thrust  (6) .  With  the  associated 
power  processing  units  (PPU)  and  interface  hardware,  the 
mass  of  this  propulsion  system  is  approximately  506  kg. 
Again,  no  radiators  are  considered,  since  the  PPUs  are 
assumed  to  provide  adequate  radiation  exchange.  With  an 
additional  10%  for  electronics  and  PPU  losses,  the  ion 
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transfer  vehicle  requires  60  kW,  making  the  power  system 
mass  approximately  900  kg.  Thus,  the  dry  mass  of  the  ion 
transfer  vehicle  is  approximately  1406  kg,  again  excluding 
fuel  tanks  and  structure . 

TVS  Mass  Calculations ■  To  calculate  the  desired 
parameters  we  must  first  estimate  the  vehicle  mass 
breakdown.  For  the  first  leg  of  the  mission,  the  initial 
mass  of  the  TVS  leaving  LEO  is  given  by 


m=m  +  m  +m  +m 

IP  P  L  V 

GEO  LEO 


CA-1) 


where 


m 


GEO 


m 


LEO 


=  Mass  of  Propellant  required  for  transfer  to  GEO 

s  Mass  of  Propellant  required  for  return  transfer 

to  LEO 

=  Mass  of  Payload 

=  Mass  of  Transfer  Vehicle  (propulsion  system, 

power  system,  fuel  tanks,  structure,  and 
electronics  package) 


After  releasing  the  payload  at  GEO,  at  the  beginning  of 
the  return  leg  to  LEO,  the  total  mass  is 


m 


m 


m  =  m 

R  P 


+  m 


(A-2) 


LEO 


The  propellant  masses  can  be  found  from  the  relation 
(11:137) 


A-3 


Av  =  I  g  1 
SP  c 


n  f — ~ — 1 

Im  -m  ] 

N.  T  p' 


(A-3) 


where  g is  the  gravitational  constant  <9.81  m/s  )  and  Av  is 
the  total  change  in  velocity  required  between  the  two 
orbits.  This  velocity  change  can  be  obtained  using  the 
Edelbaum  approximation  (10:463)  for  a  continuous  spiral 
trajectory  between  two  orbits 


/2  2 

v  +  v  -  2 v  v  cos 
1  2  12 


(A-4) 


which,  for  planar  transfers  (Ai=0) ,  simplifies  to 


Av  =  v  -  v 
1  2 


(A-5) 


However,  the  total  velocity  change  required  using  the 
constrained  radii  control  law  discussed  by  this  study  is 
nearly  1.75  times  higher  than  that  given  by  the  spiral 
transfer.  Thus,  the  total  Av  will  be  increased  by  this 
factor  in  further  calculations . 

The  velocity  of  the  transfer  vehicle  in  circular  orbit 
about  the  Earth  is  given  by  (2:34) 


r  _  + 


(A-6) 


where 


u  =  3.986012  x  10 


s  km' 


r^.  =  6378.145  km 


h  =  altitude 
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Using  the  values  of  the  altitude  at  LEO  and  GEO,  the 
velocities  are 


v  =  7725.76  m/s 

LEO 

v  =  3071 .86  m/s 

OEO 


and  the  magnitude  of  the  velocity  change  required  to 
transfer  between  LEO  and  GEO,  including  the  scale  factor  for 
the  constrained  radii  transfer,  is 


|Av| 


1.75  lv  -  v  I 

1  LEO  OEO 1 

8144.33  m/s  (4653.90  m/s) 


where  the  value  in  parenthesis  is  for  the  spiral  transfer. 
•  Solving  equation  (A-3)  for  the  mass  of  the  propellant 


1  exp  (i 1  "4"]  ] 

m  =  m 

P  T 

• 

SP  J 

(A-7) 


Rewriting  this  equation  in  terms  of  the  masses  at  the 
beginning  of  the  initial  transfer  to  GEO 

mp  =  mi  [  1  ~  expfj  L4zll  1  (A-8) 

OEO  ^  v  SP  J 


Assumi 
arc j  et 


ng  an 
TVS, 


initial  mass  of  5000  kg  for  both  the  ion  and 
this  equation  yields 


A-5 


m“rc  =  3089  kg  <2115  kg) 

GEO 

mlon  =  884  kg  (526  kg) 

GEO 

If  an  additional  10.5%  of  the  propellant  mass  is  assumed  to 
be  needed  for  the  fuel  tanks  and  structure  of  the  arcjet 
propulsion  package,  and  14.5%  for  the  ion,  the  transfer 
vehicle  dry  masses  increase  to 


m^rc  =  475  kg  +  .105  (3089  kg) 

i 

=  799  kg  (697  kg) 

m^°n  =  1406  kg  +  .145  (884  kg) 

l 

=  1534  kg  (1482  kg) 


The  difference  between  the  two  mass  fractions  arises  due  to 
the  additional  tank  structure  required  to  support  the 
pressurized  Xe . 

Rewriting  equation  (A-7)  in  terms  of  the  propellant 
required  for  the  return  trip  to  LEO  and  substituting 
equation  (A-2)  yields 
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Since  the  additional  tank  mass  required  to  support  the 
return  fuel  has  not  yet  been  included  in  the  vehicle  mass, 
equation  A-9  must  be  rewritten  as 


m 


LEO 


mv  +  x  m  exp[r^]  "  1 

1  rLEO  J  v  SP^c1' 


m 

V 

i 

r 

exp 

■ _ 

|Av| 
I  g 

k  SP3C 

I'1 

1  -  X 

exp^  J]  -  1 

where  the  parameter  x  is  equal  to  .105  for  the  arcjet  system 
and  .145  for  the  ion  system. 

Solving  equation  (A-10)  for  each  system  provides  return 
propellant  masses  of 

rt£rc  =  1557  kg  (553  kg) 

LEO 

n£on  =  340  kg  (177  kg) 

LEO 

making  the  final  total  transfer  vehicle  dry  masses 

799  kg  +  .105  (1557  kg) 

963  kg  (755  kg) 

1534  kg  +  .145  (340  kg) 

1584  kg  (1508  kg) 

The  only  unknown  masses  remaining  are  the  those  of  the 


m. 


m 


payloads.  Solving  equation  (A-l)  for  the  payload  mass 


(A-ll) 


m  —  m  —  m  —  m  —  m 

LIP  P  V 

OEO  LEO 


Substituting  for  the  arcjet  system 


m 


=  -609  kg 


(1577  kg) 


and  the  ion  system 


mlon  =  2192  kg 


(2789  kg) 


Thus,  it  is  obvious  that  the  constrained  radii  transfer, 
with  its  much  larger  Av  requirement,  cannot  be  performed  by 
the  arcjet  TVS. 

Table  (A-l)  provides  a  summary  of  the  final  breakdown 
of  mass  for  each  of  the  two  5000  kg  transfer  vehicle 
systems.  With  these  quantities  known,  the  maximum  percent 
change  in  the  vehicle  mass  over  one  orbit  can  now  be 
estimated.  This  maximum  change  should  occur  at  GEO,  after 
the  payload  has  been  released,  and  the  transfer  vehicle  is 
beginning  its  return  trip  to  LEO.  At  this  point,  the 
orbital  period  is  maximum  (24  hr)  and  the  TVS  total  mass 
reduced  to  the  mass  of  the  transfer  vehicle  and  the  fuel 
required  for  the  return  trip. 

The  mass  expelled  from  the  TVS  over  this  orbit  can  be 
calculated  knowing  the  mass  flow  rate  (constant)  and  the 
orbital  period  (TTIP)  .  The  mass  flow  rate  of  each  system  is 
given  by  (11:29) 

F 


m  = 


g  I 

SP 


(A-12) 


A-8 


where  F  is  the  total  vehicle  thrust 


Using  this  equation 


the  mass  flow  rate  of  each  TVS  is 


•  arc 

m 


_ 3.367  N _ 

(9.81  m/s2) (863  sec) 


.398  x  10  3  kg/sec 


•  von 
m 


5  (389  mN) _ 

(9.81  m/s2)  (4267  sec) 


=  .046  x  10  3  kg/ sec 


TVS  MASS  BREAKDOWN 

TRANSFER  VEHICLE  j 

ARCJET 

ION 

Transfer  Vehicle  Mass 

P<S9  kg 
<755  kg> 

1584  kg 
<  1508  kg> 

Fue  l  Ma  s  s  (Transfer  to  <3  EO) 

* 

<2115  kg  > 

884  kg 

<  52  <S  kg> 

Fuel  Mass  (Return  to  LEO) 

* 

<553  kg  > 

340  kg 
<177  kg> 

Paylaod  Mass 

* 

<1577  kg  > 

2192  kg 
<  2789  kg> 

Total  Mass 

5000  kg 

5000  kg 

Table  A-l .  Transfer  Vehicle  System  Mass  Distribution 
required  for  a  Constrained  Radii  of 
Perigee  and  Apogee  Transfer 


Note:  ()  refer  to  quantities  for  a  spiral  transfer 

*  Arcjet  transfer  vehicle  incapable  of  mission 
using  constrained  radii  control  law 


The  change  in  mass  (mass  expelled)  over  one  orbit  can 
be  obtained  using 

Am  =  m  TTIP  (A- 13) 

where  TTIP  is  the  period  of  one  orbit.  Thus,  for  each 
propulsion  system 

Am"0  =  marC  TTIP  =  (.000398  kg/sec)  (24  hr)  (3600  sec/hr) 

=  34.36  kg 

Am''or’  =  (.000046  kg/sec)  (24  hr)  (3600  sec/hr) 

=  4 .01  kg 

giving  a  final  estimate  of  the  maximum  percent  change  in 

masses  over  one  orbit  as  (spiral  transfer  only  for  arcjet 
TVS) 

^  (10096)  =  (l3Q86kg?  <100*>  =  (2.6396) 


A  v  on 

—  (100%) 

mR 


4.01  kg 
1924  kg 


(100%) 


.2196  (.24%) 


Maximum  Acceleration  Calculations ■  The  acceleration  of 
the  TVS  is  given  by 

A  =  —  (A-14) 

m 

For  comparison,  we  will  calculate  the  anticipated  maximum 
acceleration  on  the  TVS  at  both  LEO  and  GEO. 

LEO .  The  maximum  acceleration  experienced  by  the 
TVS  during  the  mission  occurs  upon  final  return  of  the  TVS 


A-l  0 


to  LEO.  At  this  point,  since  all  fuel  has  been  consumed, 
the  TVS's  mass  has  been  reduced  to  that  of  the  transfer 
vehicle  dry  weight.  Thus,  the  instant  the  last  of  the  fuel 
is  consumed,  the  acceleration  on  the  arcjet  system  (spiral 
transfer  only)  is 


arc 


-K) 


3.367  N  .  .  2 . 

=  (.00446  m/s  ) 


(755  kg) 


-  s 


=  (.45  x  10  g's) 


and  the  ion  system 


F 

m 


5 15  84^kg^~  =  •  00123  ™/s2  C-  00129  m/s2) 


=  .13  x  10"3  g’s  (.13  x  10  3  g's) 


GEO ■  The  largest  acceleration  experienced  at  GEO 
occurs  just  after  the  release  of  the  payload,  when  the  TVS 
begins  the  final  return  leg  of  the  mission.  For  the  arcjet 
propulsion  system,  this  acceleration  (spiral  only)  is 


m 


arc 


3.367  N  =  ( t Q0257  m/s2) 


(1308  kg) 


=  (.26  x  10"3g'e) 


and  the  ion  system 


_  V  Ol 


F 

m 


5  (389  mN)  _  rtri,  y  2  ,  ,  2, 

^ ^24 — ^ —  .00101  m/s  (.00115  m/s  ) 


=  .10  x  10‘3  g’s  (.12  x  10-3  g’s) 
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TVS  Transfer  Times .  For  a  continuous  thrust  propulsion 
system,  the  mass  flow  rate  is  constant  allowing  the  transfer 
time  to  be  calculated  using  the  relation 


t 


exp 

m 


(A-15) 


where  m  is  the  total  mass  expelled.  Thus,  the  transfer 

exp 

times  for  each  TVS  can  be  estimated  by  simply  dividing  the 
total  fuel  mass  expelled  on  each  leg  of  the  mission  (LEO  to 
GEO,  and  return  to  LEO)  by  the  total  mass  flow  rate.  This 
results  transfer  times  to  GEO  of 


r  m 

r  ~  arc  I 

p 

a  eo 

i 

"o' 

u 

0 

J 

.  m 

(2114 

kg) 

* 

1  hr 

,  -  ^ 
1  day 

.000398 

kg/ sec 

3600  sec 
»  * 

24  hrs 
•  * 

(61.53  days,  spiral  only) 


r  m  -x 

r  ^ i on  I 

P  1 

L  1  _  1 

OEO 

i 

— s 

0 

U 

0 

j 

.  m  J 

884  kg 

f  1  hr 

1  day  1 

.000046  kg/sec 

|  3600  sec 

24  hrs 

t.  - 

=  220.21  days 

(131.05 

days) 
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and  return  to  LEO 


ft  1 

^  LEO J 


m 


LEO 


m 


(55  3  kg) 
000398  kg/sec 


C  f 

1  hr  1 

3600  sec  24 

i.  /  V. 


day 


hrs 


=  (16.10  days,  spiral  only) 


(wT-  ( 


m 


LEO 


m 


ion 


340 

kg 

f  1  hr 

■v 

* 

1 

day 

000046 

kg/ sec 

3600  sec 

J 

24 

hrs 

j 

=  84.72  days 


(44 . 17  days) 


Summary  of  Results ■  Table  (A-2)  provides  an  outline  of 
the  propulsion  system  and  vehicle  performance  parameters, 
mass  distribution,  and  transfer  times  for  the  N  H  arc jet 

Z  4 

and  Xe  ion  transfer  vehicles.  As  a  comparison,  this  data 
includes  the  results  of  utilizing  the  constrained  radii  of 
apogee/perigee  transfer  and  the  spiral  transfer. 
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** 

— 

1 

SYSTEM 

CONSTRAINED  RADII  | 

SPIRAL 

PARAMETERS 

arcjet 

i  on 

arc  jet 

i  on 

Thruster  Perform. 

Specif ic  Impulse 

8  <53  sec 

42<5?  sec 

803  sec 

4  207  sec 

Eff  iciency 

35  % 

75  96 

35  * ** 

75  % 

Power/TSruster 

2<5.  O  kW 

10.87  kW 

2<5.  O  kW 

10.87  kW 

Thrust/Thruster 

3.  3  <57  N 

380  mN 

3. 3 <57  N 

380  mN 

#  of  Thrusters 

1 

5 

1 

5 

Power  System 

Mas  s/Pover 

15  kg/kw 

15  kg/kw 

15  kg/kW 

15  kg/kw 

Ma  ss  Distribution 

Propulsion  Sys 

25  kg 

50<5  k  g 

25  kg 

500  kg 

Power  System 

450  kg 

0OO  k  g 

450  kg 

POO  k  g 

Tanks/Structure 

* 

178  kg 

280  kg 

102  kg 

Fue  l 

4» 

1224  kg 

2<5<58  kg 

703  k  g 

Pay  load 

* 

2102  kg 

1577  kg 

2780  kg 

Tot  al 

* 

5000  kg 

5000  kg 

5000  kg 

Max  xAm/orbit 

* 

O.  21  M 

2.  <53  96 

0.24  96 

Mass  Flow  Rate 

. 000308 

.  00004(5 

.  000308 

.  00004  0 

<  km/sec ) 

Max  Acceleration 

LEO  , 

41 

. OOOl 3 

.  00045 

.  00013 

<  q  s  > 

OEO  3 

* 

. OOOIO 

.  0002  <5 

.  00012 

Transfer  Time 

1 

* 

220  days 

<S2  days 

131  days 

41 

85  days 

lO  days 

44  day  s 

Table  A-2 .  Performance  Parameters  and  Mass  Distribution  o 
Ammonia  Arcjet  and  Xenon  Ion  Transfer  Vehicles 
for  Constrained  Radii  and  Spiral  Tranfers. 


*  Arcjet  TVS  incapable  of  constrained  radii  transfer 

**  Total  mission  includes  transfer  of  payload  from  LEO 
(300  km)  to  GEO  (35,863  km)  and  return  of  "empty" 
transfer  vehicle  to  LEO. 
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Appendix  B-l:  Program  DELAMAX2 


C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*»***»»*»**»»*»»»»»*#*»*»#*»»** 


ft  ft 

*  PROGRAM  DELAMAX2  * 
«  * 
«  Capt  Gregory  Beeker  * 

*  Air  Force  Institute  of  Technology  * 
»  May  1988  * 
»  « 


ft****************************** 

The  following  program  solves  for  the  maximum  change  in  the 
magnitude  of  the  nondimens ional  semimajor  axis  for  the  contin¬ 
uous  thrust  orbit  transfer  problem  in  which  the  distance  to 
perigee  (Rp)  or  apogee  (Ra)  is  held  fixed.  Simpson's  rule 
is  used  in  subroutine  integrate  to  solve  the  integral 
equations  for  the  change  in  the  nondimensional  semimajor 
axis  (delta  a*)  and  eccentricity  (delta  e»)  given  values  of 
a,  e,  nu,  and  the  Lagrange  multiplier,  lambda.  The  secant 
method  is  used  to  find  the  value  of  lambda  which  will  drive 
the  value  of  delta-Rp  or  delta-Rp  to  zero. 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) ,NUI (361) 

COMMON  E, DELR.DELA, DELE, LAMBDA, IDIV, SIGN, ALPHAI ,NUI 

*  INPUT  INITIAL  DATA  # 

«»»»»»»»»»##»*»»«*»»»» 

WRITE  («,*)  ’Please  enter  the  initial  and  final  values’ 

WRITE  (*,*)  ’of  the  eccentricity,  e,  dear.’ 

READ  ( *  ,  » )  EO  ,  EF 

WRITE  f»,»)  'Enter  the  number  of  steps  between  the  initial’ 

WRITE  (»,»)  ’and  final  values  of  eccentricity.’ 

READ  (»,*)  INUME 

WRITE  (»,»)  ’Do  you  wish  to  keep  the  distance  to  ’ 

WRITE  (»,»)  'apogee  constant,  or  the  distance  to’ 

WRITE  (*,*)  'perigee  constant  ?’ 

WRITE  (*,*)  ’(Type:  -1.0  for  perigee;  1.0  for  apogee)’ 

REAv  (*,*)  SIGN 

WRITE  (*,»)  'Enter  the  number  of  pieces  Subroutine’ 

VRITE  (»,»)  'Integrate  is  to  dissect  the  integals  into.' 

WRITE  (#,»)  '(Must  be  an  even  number,  360  maximum)’ 

READ  (»,»)  IDIV 

WRITE  (*,*)  'Enter  \.he  initial  and  second  guess  of  the  ' 

WRITE  (»,*)  'parameter  lambda.’ 

READ  (»,*)  LAMBDAO .LAMBDA  1 
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WRITE  (*,*)  'Enter  the  maximum  number  of  iterations  to  be’ 
WRITE  (*,*)  ’performed  by  secant  method  routine.’ 

READ  (*,*)  IMAX 

WRITE  (*,*)  ’Enter  the  tolerance  of  delta-Ra/p. ’ 

READ  (*,*)  TOLR 

WRITE  (*,*)  ’What  type  of  print  out  do  you  wish?’ 

WRITE  (*,*)  ’(0  for  document  data,  1  for  plot  data)’ 

READ  (»,*)  IPNT 

WRITE  (*,*)  ’If  document  data  was  selected,  do  you  wish’ 

WRITE  (*,*)  ’to  print  the  thrust  angle  around  the  orbit?’ 

WRITE  (*,*)  ’ (0  -  no,  1  -  yes) ’ 

READ  (*,*)  ITANG 

*  Print  Initial  Data  * 
a********************* 

IF  KIPNT.EQ.l) .AND. (SIGN. LE.O.))  THEN 

OPEN  (UNIT2 1 1 ,  STATUS  2  ’NEW’, FILE  =  ’APLAM.DAT’ ) 

OPEN  (UNIT2 12 ,  STATUS  2  ’NEW', FILE  =  ’APDELA.DAT’ ) 

OPEN  (UNIT2 13 ,  STATUS  2  ’NEW’, FILE  =  ’APDELE.DAT’ ) 

END  IF 

IF  ( (IPNT. EQ. 1) .AND. (SIGN. GT.O. ) )  THEN 

OPEN  (UNIT2 1 1 ,  STATUS  2  ’NEW’, FILE  2  ’AALAM.DAT’) 

OPEN  (UNIT2 12,  STATUS  2  ‘NEW’, FILE  2  ’AADELA.DAT’ ) 

OPEN  (UNIT2 13 ,  STATUS  2  ’NEW’, FILE  2  ’AADELE.DAT’) 

END  IF 

IF  (IPNT.EQ.O)  THEN 

OPEN  (UNIT2 10 ,  STATUS  2  ’NEW’, FILE  2  ’ AMAX2 . OUT ’ ) 

WRITE  (10,12)  IMAX 
WRITE  (10,15)  IDIV 
WRITE  (10,*) 

IF  (SIGN. GT.O)  WRITE  (10,18)  TOLR 
TF  (SIGN. LT . 0)  WRITE  (10,19)  TOLR 
WRITE  (10,*) 

WRITE  (10,25) 

END  IF 

WRITE  (6,12)  IMAX 
WRITE  (6,15)  IDIV 
WRITE  (6,*) 

IF  (SIGN. GT.O)  WRITE  (6,18)  TOLR 
IF  (SIGN. LT . 0)  WRITE  (6,19)  TOLR 
WRITE  (6,*) 

WRITE  (6,25) 

C 

12  FORMAT  (3X,’Max  *  of  iterations  to  be  performed  ’, 

+  ’by  Secant  Method  Routine  is  ’,14) 

15  FORMAT  (3X, ’Number  of  divisions  for  each  integral  is  ’,14) 

18  FORMAT  (3X, ’Apogee  Distance  Fixed  with  a  tolerance  of  , E8 . 1 ) 

19  FORMAT  (3X, ’Perigee  Distance  Fixed  with  a  tolerance  of  , E8 .  1 ) 
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25  FORMAT  (6X, ’e’ . 13X, 'Lambda' , 16X, 'Delta  a# ' , 14X, ' Delta  e*’) 


ESTEP= (EF-EO) /INUME 
E=E0 


*  BEGIN  INITIAL  LOOP  WHICH  INCREMENTS  ECCENTRICITY  # 
########*##*#*##**##****####*#*###*##**#»#**#####»#* 


L AMBD A= L AMBD AO 
CALL  INTEGRATE 

IF  (ABS(DELR) . LE.TOLR)  GO  TO  45 
DELRM1=DELR 
LAMDM1= LAMBDA 
LAMBDA-LAMBDA 1 
30  CALL  INTEGRATE 

IF  (ABS(DELR) .LE.TOLR)  GO  TO  45 

*  BEGIN  ITERATIONS  OF  SECANT  METHOD  » 
*############**#*#**##*»*##**#***#*** 


C 


C 


C 

c 

c 

c 

c 


DO  40  1=1, IMAX 

LAMDP1=LAMBDA-DELR* (LAMBDA-LAMDM1) / (DELR-DELRM1) 

DELRM1=DELR 

LAMDM1= LAMBDA 

LAMBDA=LAMDP1 

CALL  INTEGRATE 

IF  (ABS(DELR) .LE.TOLR)  GO  TO  45 
40  CONTINUE 

IF  (IPNT.EQ.O)  THEN 
WRITE  (10.*) 

WRITE  (10,42)  E.TOLR 
WRITE  (10,*) 

END  IF 

WRITE  (6 , ») 

WRITE  (6,42)  E.TOLR 
WRITE  (6,*) 

42  FORMAT  (3X,’e  =  ’ ,F8 . 4 , 3X Secant  Method  did  not  converge’, 
+  ’  within  given  tolerance  of  ’  , E 1 5 . 7 ) 

GO  TO  60 


*  PRINT  FINAL  DATA  » 

a#*#####*#*****##*#* 

45  WRITE  (6,48)  E , LAMBDA , DELA , DELE 

IF  (IPNT.EQ.O)  WRITE  (10,48)  E , LAMBDA , DELA . DELE 
IF  (IPNT.EQ.l)  THEN 
WRITE  (11,49)  E, LAMBDA 
WRITE  (12,49)  E.DELA 
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WRITE  (13,49)  E.DELE 
END  IF 

IF  ( (IPNT.EQ.O) .AND. (ITANQ.EQ. 1) )  THEN 
WRITE  (6 ,  *) 

WRITE  (10,*) 

WRITE  (6,50) 

WRITE  (10,50) 

WRITE  (6,*) 

WRITE  (10,*) 

DO  46  1=1 . (IDIV+1) 

WRITE  (6,51)  NUI(I),  ALPHAI(I) 

WRITE  (10,51)  NUI (I) ,  ALPHAI(I) 

46  CONTINUE 
WRITE  (6 , ») 

WRITE  (10,*) 

END  IF 

48  FORMAT  (3X,F6 . 4 , 3X ,E20 . 13 , 3X , E20 . 13 , 3X , E20 . 13) 

49  FORMAT  (3X.E20 . 13 ,3X,E20 . 13) 

50  FORMAT  (7X, 'NU' ,8X, 'ALPHA' ) 

51  FORMAT  (3X,F7 . 2 , 5X,F7 . 2) 

DIFF=EF-E 

IF  (ABS(DIFF) .LE. IE-10)  GO  TO  60 

E=E+ESTEP 

GO  TO  30 

60  STOP 

END 


********************************************************** 

SUBROUTINE  INTEGRATE 

********************************************************** 

SUBROUTINE  INTEGRATE 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) , NUI ( 36 1 ) 

COMMON  E.DELR.DELA, DELE, LAMBDA, IDIV, SIGN, ALPHAI ,NUI 

NU=0 . 

AINT1  =  0 . 0 
AINT2=0 . 0 
EINT=0 . 0 
ESQU= 1 . -E**2 
PI=DBLE (ACOS (-1.0)) 

WRITE  (*,»)  PI 
DELNU=2 . *PI/IDIV 
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Fl  =  2 . »E+LAMBDA* (2 . #E* ( 1 . +SIGN«E) +SIGN*ESQU) 

F2=SIGN#LAMBDA»ESQU#»2 
WRITE  (»,*)  FI ,F2 

DO  100  1=1,  (IDIV+1) 

CNU= 1 . +E»DC0S (NU) 

SNU=DSIN ( NU) 

*  CALCULATE  THRUST  VECTOR  ANGLE  » 
*##«*#######*»#####**###########* 

NUMER=E«SNU»CNU*F1 
DEN0M=CNU*#2#F 1 -F2 
WRITE  (*,*)  NUMER, DENOM 

IF  ((DABS (NUMER) .LE. 1 .E-15) . AND. (DABS (DENOM) 

+  .LE. 1 .E-15) )  THEN 

WRITE  (6,80) 

IF  (IPNT.EQ.O)  WRITE  (10,80) 

80  FORMAT  (3X,’The  denominator  argument  of  the’, 

+  ’ARCTAN  function  was  equal  to  or  near  0.’) 

WRITE  (6,85)  NUMER.DENOM.Fl ,F2,NU 
IF  (IPNT.EQ.O)  WRITE  (10,85)  NUMER.DENOM.Fl ,F2 ,NU 
85  FORMAT  (3X,’NUM=  ’ .E15.8 ,3X, ’DEN  =  \E15.8,3X, 

♦  ’FI  =  ' .E15.8.3X, *F2  =  ’.EIS.S.SX, 

+  'NU  =  ’ , E15 . 8  ,  ’  Rad’) 

ALPHA=0 . 0 
GO  TO  88 
END  IF 

ALPHA=DATAN2 (NUMER , DENOM) 

88  ALPHAI (I) =180. *ALPHA/PI 
NUI (I) =180. *NU/PI 

IF  (ALPHAI (I) .LT.O.)  ALPHAI (I) =360. +ALPHAI (I) 
CALPHA=DCOS (ALPHA) 

SALPHA=DSIN( ALPHA) 


»  CALCULATE  INTEGRALS  OF  DELTA  A  AND  DELTA  E  * 
•a*##*#######***######*######*###**#*########* 


C 


C 

100 


C0NST=2 . 0 

IF  ( (1/2*2) . EQ . I )  C0NST=4 . 0 

IF  ((I.EQ. 1) .OR. (I.EQ. (IDIV+1)))  C0NST=1.0 

AINT1=AINT1+C0NST»SNU/CNU»*2»SALPHA 

AINT2=AINT2+C0NST«CALPHA/CNU 

EINT=EINT+C0NST»CALPHA/CNU«*3 

NU=NU+DELNU 

CONTINUE 
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AINTl=DELNU/3 . *AINT1 
AINT2=DELNU/3.#AINT2 
EINT=DELNU/3 . »EINT 

«  CALCULATE  VALUES  OF  DELTA-A,  DELTA-E,  AND  DELTA-Ra/p  » 
a#**#*##***#*#####*#*########**##**#####**##*#######**## 


DELA=2 . *ESQU* (E#AINT1+AINT2) 
DELE=ESQU/2 . /E*DELA-ESQU*»3/E*EINT 
DELR=DELA» ( 1 . +SIGN*E) +SIGN*DELE 

RETURN 

END 
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Appendix  B-2:  Program  DELEMAX2 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 
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»  • 

»  PROGRAM  DELEMAX2  * 

*  * 

»  Capt  Gregory  Beeker  » 

*  Air  Force  Institute  of  Technology  * 

*  May  1988  * 

«  » 


**»»««•««««««•«««««»««««««««««« 

The  following  program  solves  for  the  maximum  change  in  the 
magnitude  of  the  nondimensional  eccentricity  for  the  contin¬ 
uous  thrust  orbit  transfer  problem  in  which  the  distance  to 
perigee  (Rp)  or  apogee  (Ra)  is  held  fixed.  Simpson’s  rule 
is  used  in  subroutine  integrate  to  solve  the  integral 
equations  for  the  change  in  the  nondimensional  semimajor 
axis  (delta  a*)  and  eccentricity  (delta  e*)  given  values  of 
a,  e,  nu,  and  the  nondimensional  Lagrange  multiplier,  lambda*. 
The  secant  method  is  used  to  find  the  value  of  lambda*  which 
will  drive  the  value  of  delta-Rp  or  delta-Rp  to  zero. 


IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) ,NUI (361) 

COMMON  E , DELR , DELA , DELE , LAMBDA , IDI V , SIGN , ALPHAI , NUI 

*  INPUT  INITIAL  DATA  * 

»»*»»»»»»»»*«*»»»*«*»» 

WRITE  (*,«)  ’Please  enter  the  initial  and  final  values’ 
WRITE  (»,*)  'of  the  eccentricity,  e,  dear.’ 

READ  («, »)  EO.EF 

WRITE  (*,«)  'Enter  the  number  of  steps  between  the  initial’ 

WRITE  (»,*)  'and  final  values  of  eccentricity.’ 

READ  (»,*)  INUME 

WRITE  (*,*)  ’Do  you  wish  to  keep  the  distance  to  ’ 

WRITE  (*,»)  'apogee  constant,  or  the  distance  to’ 

WRITE  (».*)  'perigee  constant  ?’ 

WRITE  (»,»)  ’(Type:  -1.0  for  perigee;  1.0  for  apogee)’ 

READ  (»,*)  SIGN 

WRITE  (»,*)  'Enter  the  number  of  pieces  Subroutine’ 

WRITE  (»,*)  ’Integrate  is  to  dissect  the  integrals  into.’ 
WRITE  (»,*)  ’(Must  be  an  even  number,  360  maximum)’ 

READ  (»,*)  IDIV 

WRITE  (*, »)  'Enter  the  initial  and  second  guess  of  the  ’ 
WRITE  (», «)  'parameter  lambda*  (Lagrange  multiplier,’ 

WRITE  («, «)  'lambda,  times  the  semimajor  axis,  a).’ 

READ  (*, »)  LAMBDAO, LAMBDA 1 
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WRITE  (*,*)  ’Enter  the  maximum  number  of  iterations  to  be’ 
WRITE  (*,*)  ’performed  by  secant  method  routine.’ 

REAP  (*,*)  IMAX 

•  WRITE  (*,*)  ’Enter  the  tolerance  of  delta-Ra/p. ’ 

READ  (*,*)  TOLR 

WRITE  (*,*)  ’What  type  of  print  out  do  you  wish?’ 

WRITE  (*,*)  ’(0  for  document  data,  1  for  plot  data)’ 

READ  (*,*)  IPNT 

WRITE  (*,*)  ’If  document  data  was  selected,  do  you  wish’ 

#  WRITE  (*,*)  ’to  print  the  thrust  angle  around  the  orbit?’ 

WRITE  (*,*)  ’ (0  -  no,  1  -  yes) ’ 

READ  (*,*)  ITANG 

*  Print  Initial  Data  * 

********************** 


IF  ( (IPNT. EQ. 1) .AND. (SIGN. LE.O.) )  THEN 
OPEN  (UNIT= 1 1 ,  STATUS  =  ’NEW’, FILE  =  ’PLAM.DAT’ ) 

OPEN  (UNIT= 12 ,  STATUS  =  ’NEW’, FILE  =  ’PDELA.DAT’ ) 

OPEN  (UNIT* 13 ,  STATUS  =  ’NEW’ .FILE  =  ’PDELE.DAT’) 

END  IF 
C 

IF  ( (IPNT.EQ. I) .AND. (SIGN.GT.O. ) )  THEN 
OPEN  (UNIT= 1 1 ,  STATUS  =  ’NEW* .FILE  =  ’ALAM.DAT’) 

OPEN  (UNIT= 12 ,  STATUS  =  ’NEW’ .FILE  =  ’ADELA.DAT’) 

OPEN  (UNIT= 13 ,  STATUS  =  ’NEW’, FILE  =  ’ADELE.DAT’) 

END  IF 
C 

IF  (IPNT.EQ.O)  THEN 

OPEN  (UNIT= 10 ,  STATUS  =  ’NEW’ .FILE  =  ’ EMAX2 . OUT ’ ) 

WRITE  (10,12)  IMAX 
WRITE  (10,15)  IDIV 
WRITE  (10,*) 

IF  (SIGN.GT.O)  WRITE  (10,18)  TOLR 
IF  (SIGN. LT . 0)  WRITE  (10,19)  TOLR 
WRITE  (10,*) 

WRITE  (10,25) 

END  IF 
C 

WRITE  (8,12)  IMAX 
WRITE  (8,15)  IDIV 
WRITE  (6 ,  *) 

IF  (SIGN.GT.O)  WRITE  (6,18)  TOLR 
IF  (SIGN. LT. 0)  WRITE  (6,19)  TOLR 
WRITE  (6,*) 

WRITE  (6,25) 

C 

12  FORMAT  (3X,’Max  »  of  iterations  to  be  performed  ’, 

♦  ’by  Secant  Method  Routine  is  ’,14) 

15  FORMAT  (3X, ’Number  of  divisions  for  each  integral  is  ’,14) 

18  FORMAT  (3X, ’Apogee  Distance  Fixed  with  a  tolerarce  of  ’  , E8 . 1 ) 

19  FORMAT  (3X, ’Perigee  Distance  Fixed  with  a  tolerance  of  ’  , E8 . 1 ) 
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25  FORMAT  (6X, ’e’ , 13X, ’Lambda*’ , 16X, ’Delta  a* ’ , 14X, ’Delta  e*’) 


ESTEP= (EF-EO) /INUME 
E=E0 


*  BEGIN  INITIAL  LOOP  WHICH  INCREMENTS  ECCENTRICITY  * 
••a**#*##****#*#*****####*****#**#**##*##*****##*#** 

LAMBDA=LAMBDAO 
CALL  INTEGRATE 

IF  (ABS(DELR) .LE.TOLR)  GO  TO  45 
DELRM1=DELR 
LAMDM1=LAMBDA 
LAMBDA^ LAMBDA 1 
30  CALL  INTEGRATE 

IF  (ABS(DELR) .LE.TOLR)  GO  TO  45 

*  BEGIN  ITERATIONS  OF  SECANT  METHOD  * 
a#***###*****#*#*****######**#**#*#*# 

DO  40  1=1, I MAX 

LAMDP 1 =LAMBDA-DELR* (LAMBDA-LAMDM1 ) / (DELR-DELRM1) 

DELRM1=DELR 

LAMDM1=LAMBDA 

LAMBDA* LAMDP 1 

CALL  INTEGRATE 

IF  (ABS(DELR) .LE.TOLR)  GO  TO  45 
40  CONTINUE 

IF  (IPNT.EQ.O)  THEN 
WRITE  (10, ») 

WRITE  (10,42)  E.TOLR 
WRITE  (10,*) 

END  IF 

WRITE  (6.*) 

WRITE  (6,42)  E.TOLR 
WRITE  (6,*) 

42  FORMAT  (3X, 'e  =  ’ ,F8 . 4 ,3X, ’ Secant  Method  did  not  converge’, 
+  ’  within  given  tolerance  of  ’ ,E15.7) 

GO  TO  60 


»  PRINT  FINAL  DATA  * 

»*•»»»»«****»«*»«»*« 

45  WRITE  (6,48)  E, LAMBDA ,DELA, DELE 

IF  (IPNT.EQ.O)  WRITE  (10,48)  E , LAMBDA , DELA , DELE 
IF  (IPNT.EQ.l)  THEN 
WRITE  (11,49)  E, LAMBDA 
WRITE  (12,49)  E.DELA 
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c 


c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


WHITE  (13,49)  E , DELE 
END  IF 

IF  ( (IFNT.EQ.O) .AND. (ITANG.EQ. 1) )  THEN 
WRITE  (6,*) 

WRITE  (10,#) 

WRITE  (6,50) 

WRITE  (10,50) 

WRITE  (8,#) 

WRITE  (10,#) 

DO  46  1=1 , (IDIV+1) 

WRITE  (6,51)  NUI (I) ,  ALPHAI(I) 

WRITE  (10,51)  NUI(I),  ALPHAI(I) 

46  CONTINUE 
WRITE  (6,#) 

WRITE  (10,#) 

ENDIF 

48  FORMAT  (3X.F6 . 4 , 3X.E20 . 13 , 3X , E20 . 13 , 3X, E20 . 13) 

49  FORMAT  (3X.E20 . 13 , 3X , E20 . 13) 

50  FORMAT  (7X, ’NU’ ,8X, ’ALPHA’ ) 

51  FORMAT  (3X.F7 . 2 , 5X,F7 . 2) 

DIFF=EF-E 

IF  (ABS(DIFF) .LE. IE-10)  GO  TO  60 

E=E+ESTEP 

GO  TO  30 

60  STOP 

END 


#*##*################*#######*##»########***#*###*##*##### 

SUBROUTINE  INTEGRATE 

#»#########«#*######*#########*##**##**#*###»*############ 

SUBROUTINE  INTEGRATE 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) ,NUI (361) 

COMMON  E.DELR.DELA, DELE, LAMBDA, IDIV, SIGN, ALPHAI ,NUI 


NU=0 . 

AINT1=0 . 0 
AINT2=0 . 0 
EINT=0 . 0 
ESQU= 1 . -E»»2 
PI=DBLE ( ACOS (-1.0)) 
WRITE  (*,»)  PI 
DELNU=2 ,»PI/IDIV 


B-2.4 


ooooca  ca  o  o  o  o  o  o  no 


F1=ESQU+LAMBDA* (2 . *E* ( 1 . +SIGN*E) +SIGN»ESQU) 

F2= (1 . +SIGN*LAMBDA) *ESQU*»2 
WRITE  (*,*)  F1.F2 

DO  100  1=1,  (IDIV+1) 

CNU= 1 . +E*DC0S ( NU) 

SNU=DSIN(NU) 

*  CALCULATE  THRUST  VECTOR  ANGLE  * 

if*###*#**####*#*#**#####****#***# 

NUMER=E*SNU*CNU*F1 
DEN0M=CNU**2»F1-F2 
WRITE  (*,*)  NUMER.DENOM 

IF  ( (DABS (NUMER) .LE. 1 . E- 1 5) .AND. (DABS (DENOM) 

+  . LE. l.E-15) )  THEN 

WRITE  (6,80) 

IF  (IPNT.EQ.O)  WRITE  (10,80) 

FORMAT  (3X,’The  denominator  argument  of  the', 

+  'ARCTAN  function  was  equal  to  or  near  0.’) 

WRITE  (6,85)  NUMER, DENOM, F1,F2,NU 
IF  (IPNT.EQ.O)  WRITE  (10,85)  NUMER, DENOM, FI ,F2 ,NU 
FORMAT  ( 3X , ' NUM  =  ' ,E15 . 8 ,3X, ’ DEN  =  '.E15.8.3X, 

♦  'FI  =  ' .E15.8.3X, *F2  =  ',E15.8,3X, 

+  'NU  =  ' , E15 . 8 , '  Rad') 

ALPHA=0 . 0 
GO  TO  88 
END  IF 

ALPHA=DATAN2 (NUMER , DENOM) 

ALPHAI ( I ) = 180 . * ALPHA/PI 
NUI (I) =180. »NU/PI 

IF  (ALPHAI (I) .LT.0.)  ALPHAI (I) =360 .+ ALPHAI ( I) 
CALPHA=DC0S (ALPHA) 

SALPHA=DSIN (ALPHA) 


*  CALCULATE  INTEGRALS  OF  DELTA  A  AND  DELTA  E  » 
****#####*##**#»####»*##»#####»###*#***#»**»** 

C0NST=2 . 0 

IF  ((1/2*2) .EQ.I)  C0NST=4 . 0 

IF  ((I. EQ.I). OR. (I. EQ. (IDIV+1)))  C0NST=1.0 

AINT1=AINT1+C0NST*SNU/CNU**2*SALPHA 

AINT2=AINT2+C0NST*CALPHA/CNU 

EINT=EINT+C0NST*CALPHA/CNU**3 

NU=NU+DELNU 


CONTINUE 
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AINT1=DELNU/3.*AINT1 
AINT2=DELNU/3 . *AINT2 
EINT=DELNU/3.*EINT 


*  CALCULATE  VALUES  OF  DELTA-A,  DELTA-E,  AND  DELTA-Ra/p  * 
*##*»#**#*#*###*»**#*###*#*##**#*#*#*###*##**####*#*#»»# 


DELA=2 . »ESQU* (E*AINT1+AINT2) 
DELE=ESQU/2 . /E«DELA-ESQU»*3/E»EINT 
DELR=DELA» ( 1 . +SIGN*E) +SIGN»DELE 

RETURN 

END 
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Appendix  B-3:  Program  DELAMAXUC 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 


a##*##**####**#######*###***### 


»  * 

«  PROGRAM  DELAMAXUC  » 

#  « 

*  Capt  Gregory  Beeker  * 

*  Air  Force  Institute  of  Technology  # 

*  August  1988  » 

*  * 


*#»»****###*##»*###*#»*»#*#*##* 

The  following  program  solves  for  the  maximum  change  in  the 
magnitude  of  the  nondimensional  semimajor  axis  for  the  contin¬ 
uous  thrust  orbit  transfer  problem.  Simpson’s  rule  is  used 
in  subroutine  integrate  to  solve  the  integral  equations  for 
the  change  in  nondimensional  semimajor  axis  (delta  a*)  and 
eccentricity  (delta  e»)  given  values  of  a,  e,  and  nu. 


IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) ,NUI (361) 

COMMON  E.DEI  RA.DELRP.DELA, DELE, IDIV, ALPHAI ,NUI 

*  INPUT  INITIAL  DATA  * 

*#*##***»#**«#*#*#*#** 

WRITE  (#,»)  ’Please  enter  the  initial  and  final  values’ 
WRITE  (*,#)  'of  the  eccentricity,  e.’ 

READ  (#,»)  EO.EF 

WRITE  (#,»)  ’Enter  the  number  of  steps  between  the  initial’ 
WRITE  (»,»)  'and  final  values  of  eccentricity.’ 

READ  (»,*)  INUME 

WRITE  (*,»)  'Enter  the  number  of  pieces  Subroutine’ 

WRITE  (»,*)  'Integrate  is  to  dissect  the  integrals  into.’ 
WRITE  («,»)  ’(Must  be  an  even  number,  360  maximum)' 

READ  (#,#)  IDIV 

WRITE  (*,*)  'What  type  of  print  out  do  you  wish?’ 

WRITE  (*,*)  ’(0  for  document  data,  1  for  plot  data)’ 

READ  (»,#)  IPNT 

WRITE  (»,»)  'If  document  data  was  selected,  do  you  wish’ 
WRITE  (*,»)  'to  print  the  thrust  angle  around  the  orbit?’ 
WRITE  (*,*)  ’(0  -  no,  1  -  yes)’ 

READ  (*,*)  ITANG 

*  Print  Initial  Data  » 

it******##############* 

IF  (IPNT.EQ.l)  THEN 
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C 


C 


c 


c 


c 

c 


OPEN  (UNIT= 12 ,  STATUS 
OPEN  (UNIT* 13,  STATUS 
OPEN  (UNIT* 14,  STATUS 
OPEN  (UNIT= 15 ,  STATUS 
END  IF 


’NEW* .FILE  = 
’NEW* .FILE  = 
’NEW’ .FILE  = 
’NEW* .FILE  = 


’UCADELA.DAT’) 
’UCADELE.DAT') 
’UCADELRA.DAT’ ) 
*  UCADELRP . DAT ’ ) 


IF  (IPNT.EQ.O)  THEN 

OPEN  (UNIT= 10 ,  STATUS  =  ’NEW’ .FILE  =  ’ UCAMAXF . OUT ’ ) 
WRITE  (10,15)  IDIV 
WRITE  (10.*) 

WRITE  (10,25) 

END  IF 

WRITE  (0,15)  IDIV 
WRITE  (6,*) 

WRITE  (6.25) 


15  FORMAT  (3X, ’Number  of  divisions  for  each  integral  is  ’,14) 

25  FORMAT  (6X, ’e 13X, ’Delta  a* ’ , 14X, ’Delta  e» ’, 14X, ’Delta  ra» ’ , 
+  15X, ’Delta  rp*’) 


ESTEP= (EF-EO) /INUME 
E=E0 


30  CALL  INTEGRATE 


*  PRINT  FINAL  DATA  * 

******************** 

45  WRITE  (6,48)  E . DELA .DELE , DELRA , DELRP 

IF  (IPNT.EQ.O)  WRITE  (10,48)  E , DELA , DELE , DELRA, DELRP 

•  IF  (IPNT.EQ. 1)  THEN 

WRITE  (12,49)  E.DELA 
WRITE  (13,49)  E.DELE 
WRITE  (14,49)  E, DELRA 
WRITE  (15,49)  E, DELRP 
ENDIF 

•  C 

IF  ( (IPNT.EQ.O) .AND. (ITANG.EQ. 1) )  THEN 
WRITE  (6,*) 

WRITE  (10.*) 

WRITE  (6,50) 

WRITE  (10,50) 

•  WRITE  (6,^) 

WRITE  (10,*) 

DO  46  1=1 , (IDIV+1) 

WRITE  (6,51)  NUI ( I ) ,  ALPHAI(I) 

WRITE  (10,51)  NUI(I),  ALPHAI(I) 
m  46  CONTINUE 

•  WRITE  (6 , ») 
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WRITE  (10,*) 

ENDIF 

C 

48  FORMAT  (3X.F6 . 4 , 3X.E20 . 13 . 3X , E20 . 13 , 3X.E20 . 13 , 3X.E20 . 13) 

49  FORMAT  (3X.E20 . 13 , 3X.E20 . 13) 

50  FORMAT  (7X. ’NO* ,8X, ’ALPHA’ ) 

51  FORMAT  (3X,F7 . 2 , 5X, F7 . 2) 

C 

DIFF=EF-E 

IF  (ABS(DIFF) .LE. IE-10)  GO  TO  60 
C 

E=E+ESTEP 
GO  TO  30 
C 

60  STOP 
END 
C 
C 
C 
C 
C 

C  »###***###»*###*#*#*##»*###*«»*»##**»**#*#*###**#*#####*## 

C  SUBROUTINE  INTEGRATE 

C  ##*#*#*#*#*#######**#####**#*#*##***###*#*#*#«#»#*#*#*##»* 

c 

SUBROUTINE  INTEGRATE 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAK361)  ,NUI(361) 

COMMON  E , DELRA , DELRP . DELA , DELE . IDI V , ALPHA I , NUI 
C 

NU=0 . 

AINT1=0.0 
AINT2=0 . 0 
EINT=0 . 0 
ESQU=1.-E**2 
PI=DBLE(AC0S(-1 .0) ) 

C  WRITE  (*,*)  PI 

DELNU=2.*PI/IDIV 
DO  100  1=1,  (IDIV+1) 

CNU= 1 . +E*DCOS (NU) 

SNU=DSIN (NU) 


*  CALCULATE  THRUST  VECTOR  ANGLE  * 

•  a#*#*##*#*#####*#*)**##**###*#**# 


NUMER=E*SNU 

DENOM=CNU 

C  WRITE  (*,*)  NUMER.DENOM 

C 

IF  ( (DABS (NUMER) .LE. 1 .E-15) . AND. (DABS (DENOM) 
+  .LE. 1 .E-15) )  THEN 


i 
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WHITE  (8,80) 

IF  (IPNT.EQ.O)  WRITE  (10,80) 

80  FORMAT  (3X,’The  denominator  argument  of  the', 

+  ’ ARCTAN  function  was  equal  to  or  near  0.’) 

WRITE  (8,85)  NUMER.DENOM.HU 
IF  (IPNT.EQ.O)  WRITE  (10,85)  NUMER , DEN0M, NU 
85  FORMAT  (3X,’NUM=  ’ ,E15 . 8 , 3X , ’ DEN  5  \E15.8,3X, 

+  ’NU  =  ’  ,E15 . 8 , 1  Rad1) 

ALPHA=0 . 0 
GO  TO  88 
ENDIF 
C 

ALPHA5 DAT AN2 (NUMER .DENOM) 

88  ALPHAI (I ) 5 180 . *ALPHA/PI 
NUI (I) =180. *NU/PI 

IF  (ALPHAI (I) .LT.O.)  ALPHAI (I) =360 . +ALPHAT (I) 
CALPHA=DCOS (ALPHA) 

SALPHA=DSIN (ALPHA) 


*  CALCULATE  INTEGRALS  OF  DELTA  A  AND  DELTA  E  * 
********************************************** 

C0NST=2 . 0 

IF  ((1/2*2) . EQ . I )  C0NST=4 . 0 

IF  ( (I .EQ. 1) .OR. (I .EQ. (IDIV+1) ) )  C0NST=1.0 

AINT1=AINT1+C0NST*SNU/CNU**2*SALPHA 

AINT2=AINT2+C0NST*CALPHA/CNU 

EINT=EINT*CONST*CALPHA/CNU#*3 

NU=NU+DELNU 

100  CONTINUE 


AINTl=DELNU/3 . *AINT1 
AINT2=DELNU/3 . *AINT2 
EINT=DELNU/3 . *EINT 

*  CALCULATE  VALUES  OF  DELTA-A,  DELTA-E,  AND  DELTA-Ra/p  * 
**•»**«**«»****»•«**»**»»«*•»«****»»*******»*******»**** 


DELA=2 . *ESQU* (E*AINT1+AINT2) 
DELE=ESQU/2 . /E*DELA-ESQU**3/E*EINT 
DELRA=DELA» ( 1 . +E) +DELE 
DELRP=DELA* ( 1 . -E) -DELE 

RETURN 

END 
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Appendix  B-4:  Program  DELEMAXUC 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


a###*##############*########*** 


«  « 

«  PROGRAM  DELEMAXUC  « 

«  • 

*  Capt  Gregory  Beeker  * 

*  Air  Force  Institute  of  Technology  » 

*  August  1988  • 

*  » 


##*#*#**##«*######»»**##»###*## 

The  following  program  solves  for  the  maximum  change  in  the 
magnitude  of  the  nondimensional  eccentricity  for  the  contin¬ 
uous  thrust  orbit  transfer  problem.  Simpson’s  rule  is  used 
in  subroutine  integrate  to  solve  the  integral  equations  for 
the  change  in  nondimensional  semimajor  axis  (delta  a»)  and 
eccentricity  (delta  e«)  given  values  of  a,  e,  and  nu. 


IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) ,NUI (361) 

COMMON  E, DELRA.DELRP.DELA, DELE, I DIV, ALPHAI .NUI 

»  INPUT  INITIAL  DATA  » 

#*»*###*###*#*»***##** 

WRITE  (*,*)  ’Please  enter  the  initial  and  final  values’ 
WRITE  (*,*)  'of  the  eccentricity,  e,  dear.’ 

READ  (*,*)  EO.EF 

WRITE  («,*)  'Enter  the  number  of  steps  between  the  initial’ 
WRITE  (»,*)  ’and  final  values  of  eccentricity.’ 

READ  (*,*)  INUME 

WRITE  (»,»)  ’Enter  the  number  of  pieces  Subroutine’ 

WRITE  (*,*)  ’Integrate  is  to  dissect  the  integrals  into.’ 
WRITE  (»,*}  ’(Must  be  an  even  number,  360  maximum)’ 

READ  (*,*)  IDIV 

WRITE  («,*)  ’What  type  of  print  out  do  you  wish?’ 

WRITE  (»,»)  ’(0  for  document  data,  1  for  plot  data)’ 

READ  (»,*)  IPNT 

WRITE  (»,#)  'If  document  data  was  selected,  do  you  wish’ 
WRITE  (»,»)  'to  print  the  thrust  angle  around  the  orbit?’ 
WRITE  (»,*)  ’ (0  -  no,  1  -  yes)  ’ 

READ  (*,*)  ITANG 

»  Print  Initial  Data  * 

•a#######***##*****### 

IF  (IPNT. EG. 1)  THEN 
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c 


C 


C 


C 


OPEN  (UNIT=12 ,  STATUS 
OPEN  (UNIT= 13 ,  STATUS 
OPEN  (UNIT= 14 ,  STATUS 
OPEN  (UNIT=15,  STATUS 
END  IF 


'NEW' .FILE  = 
’NEW* .FILE  = 
’NEW’ .FILE  = 
’NEW’ .FILE  = 


•UCEDELA.DAT’) 
'UCEDELE.DAT' ) 
'UCEDELRA.DAT') 
'UCEDELRP.DAT') 


IF  (IPNT.EQ.O)  THEN 

OPEN  (UNIT= 10 ,  STATUS  =  ’NEW’ .FILE  =  ’ UCEMAXF . OUT ’ ) 
WRITE  (10,15)  IDI7 
WRITE  (10, *) 

WRITE  (10,25) 

ENDIF 


WRITE  (6,15)  IDIV 
WRITE  (6,») 

WRITE  (6,25) 


15  FORMAT  (3X, 'Number  of  divisions  for  each  integral  is  ’,14) 

25  FORMAT  (6X, ’ e ’.13X, 'Delta  a* 14X, 'Delta  e* ', 14X, 'Delta  ra*’, 
+  15X, 'Delta  rp*') 

ESTEP= (EF-EO) /INUME 
E=E0 


30  CALL  INTEGRATE 


*  PRINT  FINAL  DATA  * 

»###*»#*###**#*#*### 

45  WRITE  (6.48)  E , DELA , DELE , DELRA , DELRP 

IF  (IPNT.EQ.O)  WRITE  (10,48)  E , DELA , DELE , DELRA , DELRP 
IF  (IPNT.EQ. 1)  THEN 
WRITE  (12,49)  E.DELA 
WRITE  (13,49)  E.DELE 
WRITE  (14,49)  E, DELRA 
WRITE  (15,49)  E, DELRP 
ENDIF 

IF  ((IPNT.EQ.O) .AND. (ITANG.EQ. D)  THEN 
WRITE  (6,*) 

WRITE  (10,*) 

WRITE  (6,50) 

WRITE  (10,50) 

WRITE  (6 , *) 

WRITE  (10.*) 

DO  46  1=1 , (IDIV+ 1 ) 

WRITE  (6,51)  NUI(I),  ALPHAI(I) 

WRITE  (10,51)  NUI(I),  ALPHAI(I) 

46  CONTINUE 
WRITE  (6 , ») 
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WHITE  (10,*) 

ENDIF 

C 

48  FORMAT  (3X, F6 . 4 , 3X, E20 . 13 , 3X. E20 . 13 , 3X.E20 . 13 , 3X.E20 . 13) 

49  FORMAT  (3X.E20. 13 ,3X,E20. 13) 

50  FORMAT  (7X, 'MU' ,8X, 'ALPHA' ) 

51  FORMAT  (3X.F7.2.5X.F7.2) 

C 

DIFF=EF-E 

IF  (ABS(DIFF) .LE. IE-10)  GO  TO  60 
C 

E=E+ESTEP 
GO  TO  30 
C 

60  STOP 
END 
C 
C 
C 
C 
C 

C  a*####*##*##**##***#####****#*###*###**#***#####*###***### 

C  SUBROUTINE  INTEGRATE 

C  #######*#*##*#*#*#*###*»######*#»»##*###*######**##**#*##* 

c 

SUBROUTINE  INTEGRATE 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  ALPHAI (361) ,NUI (361) 

COMMON  E , DELRA , DELRP , DELA , DELE , IDI V , ALPHAI , NUI 
C 

NU=0 . 

AINT1=0 . 0 
AINT2=0 . 0 
EINT=0.0 
ESQU=1 . -E**2 
PI=DBLE (ACOS (-1.0)) 

C  WRITE  (*,*)  PI 

DELNU=2.*PI/IDIV 
DO  100  1=1,  (IDIV+1) 

CNU= 1 . +E*DCOS (NU) 

SNU=DSIN(NU) 


*  CALCULATE  THRUST  VECTOR  ANGLE  * 
*#**#»***##*****#*******#**####*» 

NUMER=E*SNU*CNU 
DENOM=CNU**2-ESQU 
C  WRITE  (»,*)  NUMER , DENOM 

C 

IF  ( (DABS (NUMER) .LE. 1 .E-15) . AND. (DABS (DENOM) 
♦  .LE.l.E-15))  THEN 
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WHITE  (6,80) 

IF  (IPNT.EQ.O)  WRITE  (10,80) 

80  FORMAT  (3X, ’The  denominator  argument  of  the', 

+  'ARCTAN  function  was  equal  to  or  near  0.') 

WRITE  (6,85)  NUMER , DENOM , NU 
IF  (IPNT.EQ.O)  WRITE  (10,85)  NUMER, DENOM.NU 
85  FORMAT  (3X,’NUM=  * , E15 . 8 , 3X, 'DEN  =  ’,E15.8,3X, 

♦  ’NU  *  ' .E15.8,  ’  Rad’) 

ALPHA=0 . 0 
GO  TO  88 
ENDIF 

ALPHA =D AT AN2 (NUMER , DENOM) 

88  ALPHAI (I) =180 . * ALPHA/PI 
NUI (I) =180. *NU/PI 

IF  (ALPHAI (I) .LT.O.)  ALPHAI (I ) =360 . + ALPHAI (I) 
CALPHA=DCOS (ALPHA) 

SALPHA=DSIN (ALPHA) 


*  CALCULATE  INTEGRALS  OF  DELTA  A  AND  DELTA  E  » 
a**###***#*##*********###***##**#########*#### 

C0NST=2.0 

IF  ((1/2*2) .EQ. I)  CONST=4 . 0 

IF  ((I.EQ.l) .OR. (I.EQ. (IUIV+1)))  CCNST=1.0 

AINT1=AINT1+C0NST*SNU/CNU**2»SALPHA 

AINT2=AINT2+CONST*CALPHA/CNU 

EINT=EINT+CONST*CALPHA/CNU**3 

NU=NU+DELNU 

100  CONTINUE 


AINT1=DELNU/3.*AINT1 
AINT2=DELNU/3.»AINT2 
EINT=DELNU/3 . *EINT 

»  CALCULATE  VALUES  OF  DELTA-A,  DELTA-E,  AND  DELTA- Ra/p  » 
»«*»«»»««*«*««»****»»*»*»»****»********«**************** 

DELA=2. *ESQU* (E*AINT1+AINT2) 

DELE=ESQU/2 . /E*DELA-ESQU**3/E»EINT 
DELRA=D£LA* ( 1 . +E) +DELE 
•  DELRP=DELA* ( 1 . -E) -DELE 

C 

RETURN 

END 
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Appendix  B-5:  Program  INTERPO 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


a***#*######*####*#*#####*##*#* 
*  « 

*  PROGRAM  INTERPO  * 

•  » 

*  Capt  Gregory  Beeker  * 

*  Air  Force  Institute  of  Technology  * 

*  June  1988  * 

«  « 

a#*#######*#****###**##***##*## 

it*##**####***#####**##**#*#**#* 

*  The  following  program  utilizes  the  Newton  formula  to  * 

*  interpolate/extrapolate  from  the  given  data  sets  * 

*  utilizing  an  nth  degree  interpolating  polynomial.  * 

*  a***#####*#**#*###*#**####*#*# 


IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E(100)  .LAMBDAUOO)  .DELA(IOO)  .DELE(IOO) 

CHARACTER  LAMBDAT  *10 

CHARACTER  DELADAT  *10 

CHARACTER  DELEDAT  *10 

CHARACTER  OUTFILE  *10 


*  INPUT  INITIAL  DATA  * 

it#*#*##*#*##**###*#### 


WRITE  (*,*)  'Please  specify  the  name  of  the  Lambda*’ 

WRITE  (*,*)  'data  file  to  be  used.’ 

READ  (* , ’ (A) ' )  LAMBDAT 

WRITE  (*,*)  'Please  specify  the  name  of  the  Delta  a*’ 

WRITE  (*,*)  'data  file  to  be  used.’ 

READ  (» , ’ (A) ' )  DELADAT 

WRITE  (*,*)  'Please  specify  the  name  of  the  Delta  e*’ 

WRITE  (*,*)  'data  file  to  be  used.’ 

READ  (» , ' (A) ’ )  DELEDAT 

OPEN  (UNIT* 15 ,  STATUS  =  'OLD' .FILE  =  LAMBDAT) 

OPEN  (UNIT* 16 ,  STATUS  =  'OLD' .FILE  =  DELADAT) 

OPEN  (UNIT* 17 ,  STATUS  =  'OLD' .FILE  =  DELEDAT) 

WRITE  (*,*)  'Is  this  data  based  on  keeping  the  distance  to  ' 

WRITE  (*,*)  'apogee  constant,  or  the  distance  to’ 

WRITE  (*,*)  'perigee  constant  ?’ 

WRITE  (*,*)  ’(0  -  apogee;  1  -  perigee)’ 

READ  (*,»)  ID 

WRITE  (*,*)  'Please  specify  the  number  of  data  points  contained’ 
WRITE  (*,*)  'in  each  of  the  data  files.’ 

READ  (*,*)  IDATA 
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WHITE  (*,*)  ’Please  specify  the  min  and  max  values  of  ' 

WHITE  (*,*)  'eccentricity  for  the  range  to  be  evaluated.’ 

READ  (*,*)  EMIN,  EMAX 

WRITE  (*  *)  'Please  enter  the  number  of  values  of  eccentricity’ 
WHITE  (*,*)  'to  be  evaluated  within  the  specified  range’ 

WHITE  (*,*)  ’(excluding  the  initial  value).’ 

READ  (*,*)  IE 

WRITE  (*,*)  'Please  specify  the  degree  of  the  interpolating’ 
WHITE  (»,*)  'polynomial  to  be  used  (cannot  exceed  n-1,  where’ 
WHITE  (*,*)  ’n  is  the  number  of  data  points.)’ 

READ  (*,*)  IDEG 
REWIND  (15) 

REWIND  (16) 

REWIND  (17) 

DO  5,  1=1, IDATA 

READ  (15, *)  E ( I ) .LAMBDA (I) 

READ  (16,*)  E(I),DELA(I) 

READ  (17,*)  Ed)  , DELE (I) 

5  CONTINUE 

WRITE  (*,*)  'Please  specify  the  name  of  the  output  file.’ 

READ  (* , ’ (A) ’ )  OUTFILE 
C 

OPEN  (UNIT3 10 ,  STATUS  =  ’NEW’, FILE  =  OUTFILE) 

C 

IF  (ID.EQ.O)  THEN 
WRITE  (6,18) 

WRITE  (10,18) 

END  IF 

IF  (ID. EQ. 1)  THEN 
WRITE  (6,19) 

WRITE  (10,19) 

END  IF 

WRITE  (6,*) 

WRITE  (10,*) 

WRITE  (6,20)  IDEG 
WRITE  (10,20)  IDEG 
C 

18  FORMAT  (3X, ’Apogee  Distance  Fixed’) 

19  FORMAT  (3X, 'Perigee  Distance  Fixed’) 

20  FORMAT  (3X, ’ Interpolating/Extrapolating  Polynomials  are  of  ’, 

+  'degree  ’ , 12) 

C 

WRITE  (6 , ») 

WRITE  (10, ») 

WRITE  (6,110) 

WRITE  (10,110) 

WRITE  (6,*) 

WRITE  (10, ») 

C 

DELTAE= (EMAX-EMIN) /IE 

IEP1=IE+1 

EBAR=EMIN 
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«  EVALUATE  FUNCTION  AT  SPECIFIED  VALUES  OF  ECCENTRICITY  « 
*«*««•«»»»«»«»»«»*»»*«»•**»*«»«*«•»•*«»«»»»»*•«***»*«*»»• 

DO  100  ICOUNT3 1 , IEP1 


*  Check  location  of  data  point  « 

DO  30  I=1,IDATA 
DIFF=EBAR-E(I) 

IF  (ABS(DIFF) .LE. l.E-15)  THEN 
LAM- LAMBDA! I) 

DA=DELA (I) 

DE=DELE(I) 

GO  TO  40 
END  IF 

IF  ((I.EQ. 1) .AND. (DIFF.LT.O.))  THEN 
IMIN= 1 

IMAX=IMIN+IDEG 
GO  TO  35 
ENDIF 

IF  ((I.EQ.IDATA) .AND. (DIFF.GT.O.))  THEN 
IMAX=IDATA 
IMIN=IMAX-IDEG 
GO  TO  35 
ENDIF 
I LEFT= I 

IF  (DIFF.LT.O.)  GO  TO  31 

30  CONTINUE 

31  IDEGP1=IDEG+ 1 
IDEGL3IDEGPl/2 

I MI N= I LEFT* 1 - I DEGL 
I MAX3 I MI N+ I DEG 
IF  (IMIN.LT. 1)  THEN 
IMIN31 

I MAX3 I MI N+ I DEG 
ENDIF 

IF  (IMAX.GT. IDATA)  THEN 
IMAXdDATA 
IMIN3IMAX-IDEG 
ENDIF 


*  Evaluate  Polynomial  • 

35  CALL  INTERP(IDLG,IMIN,IMAX,EBAR,E, LAMBDA, DELA, DELE, LAM, DA. DE) 


*  PRINT  DATA  * 
*##»#»»*##***# 
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40  WRITE  (6,120)  EBAR.LAM.DA.DE 
WRITE  (10,120)  EBAR.LAM.DA.DE 


EBAR=EBAR+DELTAE 
100  CONTINUE 


110  FORMAT  ( 13X, 'e ’ , 19X, ’Lambda* 15X, ’Delta  a»’ , 15X, ’Delta  e»’) 
120  FORMAT  (3X, E20 . 13 .3X.E20 . 13 , 3X, E20 . 13 , 3X.E20 . 13) 

CLOSE  (10) 

CLOSE  (15) 

CLOSE  (16) 

CLOSE  (17) 


*»*******»*«**»»««*«**«•«**»»**» 

*  SUBROUTINE  INTERP  * 

a#*##*#####***###**###***##*###* 

*  The  following  algorithm  evaluates  an  interpolating  * 

*  nth  degree  polynomial  * 

a#**#*#*##*#**#*#***#####*#*#**# 

SUBROUTINE  INTERP ( I DEG , I MIN , I MAX , EBAR , E . LAMBDA , 

►  DELA, DELE.LAM.DA.DE) 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E ( 100) .LAMBDA (100) .DELA(IOO) .DELE(IOO) 

DIMENSION  X(100) ,D(100) ,P(100) 


IDEGP1=IDEG+ 1 
DO  200  1=1,3 

DO  150  J=IMIN, IMAX 
K=J-IMIN+1 

IF  (I.EQ.l)  D (K) =LAMBDA(J) 

IF  (I.EQ.2)  D(K) =DELA(J) 

IF  (I.EQ.3)  D (K) =DELE ( J) 

X ( K ) =  E  ( J) 

CONTINUE 
DO  155  K= 1 , IDEG 

DO  155  J= (K+ 1) , IDEGP1 
J1=IDEGP1-J+ (K+ 1 ) 

D(J1)=(D(J1)-D(J1-1))/(X(J1)-X(J1-K)) 

CONTINUE 
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Z=EBAR-X(1) 

P ( 1 )  =D  ( 1)  +D (2)  *Z 
DO  160  J=2 , IDEG 
Z=Z* (EBAR-X(J) ) 
P(J)=P(J-1)+D(J+1)*Z 
CONTINUE 


IF  (I.EQ.l) 
IF  (I.EQ.2) 
IF  (I.EQ.3) 

CONTINUE 

RETURN 

END 


LAM=P(IDEG) 
DA=P(IDEG) 
DE=P (IPEG) 
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Appendix  C-l : 


Program  TRANSMUL 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 


*»#*#*##*#*#*#*#**#»#**#*######*« 


»  * 

*  Program  TRANSMUL  * 

»  « 

*  Captain  Gregory  Beeker  * 

*  Air  Force  Institute  of  Technology  * 

*  July  1988  » 

«  * 


*#*#»#*#*##»####*###»#»***##*###* 

The  following  program  was  written  to  solve  the  long  timescale 
problem  to  determine  the  total  transfer  delta  v  and  transfer 
time  of  a  spacecraft  traveling  between  two  planar  orbits.  The 
program  takes  a  spacecraft  from  an  initial  orbit  to  a  final 
orbit  though  many  revolutions  while  constraining  one  of  the 
two  radii  (perigee  (rp)  or  apogee  (ra)).  Subroutine  Haming  is 
incorporated  to  solve  the  differential  equations  for  the  non- 
dimensional  changes  in  eccentricity  and  semimajor  axis 
(dabar/dVbar  and  de/dVbar) .  In  addition,  Program  INTERPO,  which 
provided  the  final  solutions  to  the  fast  timescale  problem,  is 
also  incorporated  as  a  subroutine  to  provide  the  values  of  delta 
a*  and  delta  e«  as  functions  of  eccentricity. 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  I SP , MU , MANOMO , MANOM , MDOTBAR , MDOT , MO 
DIMENSION  X (4) ,VAR1(4) ,VAR2(4) ,VAR3(4) ,VAR4(4) ,VAR5(4) 
COMMON  /HAM/  VBAR, Y(42 , 4) ,F (42 , 4) , ERR (42) ,N ,H,MODE 
COMMON  /RH/  RABAR(4) ,HP6AR(4) .ICASE.ICNT, MANOMO, AO, MDOTBAR, 
+  ACCELO ,PI2 


*  INITIALIZATION  OF  PARAMETERS  FOR  HAMING  « 
»*»»»»»»**»*»»»»*»»»»»*»**»*»**»«*****»**** 


N=4 

M0DE=0 

C 

C 

C  *  INITIAL  DATA  ENTRY  * 

C  »#*##*##*#*#*#*#**#**# 

C 

C  *  Primary  Body  Parameters 

C 

WRITE  (*,»)  ’Do  you  wish  to  change  the  value  of  the’ 
WRITE  (*,*)  ’gravitational  parameter  and  radius  of  the’ 
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WRITE  (»,»)  'primary  body  (currently  set  to  those  of 
WRITE  (*,«)  'the  Earth)?  (0  -  yes;  I  -  no)’ 

READ  (»,*)  IP 

IF  (IP.EQ.O)  THEN 

WRITE  (»,»)  'Enter  the  new  values  of  mu  (km3/sec2)  and  r  (km).’ 
READ  (*,»)  M0,R 
END  IF 

IF  (IP.EQ.l)  THEN 
MU=3 . 986012D5 
R=6378 . I45D0 
END  IF 

»  Orbit  Data 


WRITE  (»,*)  ’Specify 
WRITE  (*,*) 

WRITE  (#,*)  ’ 

WRITE  (»,*)  ’ 

WRITE  (*,»)  ’ 

READ  («,*)  ITYP 


the  type  of  program  run  (0  or  1).’ 


0  -  Single  Transfer’ 

1  -  Multiple  Transfers  (Includes  plot’ 
data  file  of  af/aO  vs  delta-V  total)’ 


WRITE  (»,*)  ’Enter  the  eccentricity  and  semimajor  axis  (km)’ 
WRITE  (*,#)  'of  the  initial  orbit.’ 

READ  (*,*)  E0,A0 

2  IF  (AO.LE.R)  THEN 

WRITE ( # , *)  ’Initial  orbit  intersects  the  surface  of  the’ 
WRITE (* , *)  ’primary  body.  Please  select  a  new  value.’ 

READ  (*,»)  AO 
GOTO  2 
END  IF 


WRITE  (»,«)  ’Enter  the  eccentricity  and  semimajor  axis  (km)’ 
WRITE  («,»)  ’of  the  final  orbit.’ 

READ  (»,»)  EF.AF 

3  IF  (AF.LE.R)  THEN 

WRITE(*,«)  ’Final  orbit  intersects  the  surface  of  the’ 
WRITE (» , »)  ’primary  body.  Please  select  a  new  value.' 

READ  (»,»)  AF 
GOTO  3 
END  IF 

IF  (ITYP.EQ.P  THEN 

WRITE  (»,»)  ’Specify  the  number  of  transfers  between’ 

WRITE  (»,*)  ’the  initial  and  final  orbits.’ 

READ  (»,»)  NTRANS 
END  IF 

*  Vehicle  Parameters 
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M0=2000 
MD0T=4 . OD-4 
MDOTBAR=MDOT/MO 
ISP=5000 

WRITE  (»,»)  'Do  you  wish  to  change  the  value  of  the  vehicle  mass 
WRITE  (*,*)  'and  mass  flow  rate  parameters?  (0  -  yes;  1  -  no)’ 
WRITE  (*,*) 

WRITE  (*,*)  'Current  Settings:' 

WRITE  (*.11)  MO 
WRITE  (*,12)  MDOT 
WRITE  (*,13)  MDOTBAR 
BEAD  (*,*)  IM 

IF  (IM.EQ.O)  THEN 

WRITE  (*,*)  'Enter  the  new  values  of  the  mass  flow  rate  (kg/s) 
WRITE  (*,*)  ‘and  initial  vehicle  mass  (kg).’ 

READ  (*,»)  MDOT, MO 
MDOTBAR=MDOT/MO 
WRITE  (*,13)  MDOTBAR 
ENDIF 

WRITE  (*,*)  ’Do  you  wish  to  change  the  value  of  propulsion’ 

WRITE  (*,14)  ISP 

WRITE  (» ,*)  ’ (0  -  yes;  1  -  no)  ’ 

READ (* , *) I ISP 

IF  (IISP.EQ.O)  THEN 

WRITE  (*,*)  ’Enter  the  new  value  of  Isp.’ 

READ  (*,*)  ISP 
ENDIF 

ACCEL0  =  ISP*9 . 8 1  * MDOTBAR/ 1000. 

*  Integration  Data 

WRITE  (*,*)  ’Enter  the  step  size  of  the  independent’ 

WRITE  (*,*)  ’variable,  V-bar. ’ 

READ  (*,»)  H 

WRITE  (*,*)  ’Enter  the  maximum  number  of  iterations  to  be’ 

WRITE  (*,*)  ’performed  by  Haming  on  each  half  transfer.’ 

READ  (*,»)  IMAX 

*  Output  Parameters 


IF  (ITYP.EQ.l)  THEN 
I OUT =2 

OPEN  (UNIT  =  10,  STATUS 
OPEN  (UNIT  =11,  STATUS 
OPEN  (UNIT  =  12,  STATUS 
GO  TO  8 
ENDIF 


’NEW’,  FILE  =  ’EDELVBAR’ ) 
’NEW’,  FILE  =  ’ SDELVBAR’ ) 
'NEW’,  FILE  =  ’ HDELVBAR' ) 


o  o  o  o  o  o  o  o  o  o 


c 

WRITE  (*,»)  ’Enter  output  code.’ 

WRITE  (»,*)  *  (0  for  screen,  1  for  file  and  screen)’ 

READ  (»,«)  I0UT 
C 

8  IF  (IOUT.GE. 1)  THEM 

WRITE  («,*)  'Specify  the  name  of  the  data  output  file. 
READ  (» , ’ (A) ’ )  TRANSOUT 

OPEN  (UNIT  =  9,  STATUS  =  ’NEW’,  FILE  =  TRANSOUT) 

ENDIF 

C 

IPNT=0 

IF  (IOUT.LE. 1)  THEN 

WRITE  («,«)  ’Enter  the  number  of  steps  between  each’ 
WRITE  (*,*)  'data  printout.’ 

READ  (»,«)  IPNT 
ENDIF 
C 
C 
C 

C  *  INITIALIZE  ORBIT  PARAMETERS  ( NOND I MENS I ONAL )  » 

C  a#*#*#**##**#*##****#*#*#*####**#****#**##**##*## 

c 

C  *  State  Vector 

C 

Y2F=AF/A0 

Y2FCK=Y2F 

C 

IF  (ITYP.EQ. 1)  THEN 
DELAF  = (AF-AO) /NTRANS 
Y2F=  1 . DO 
ENDIF 
C 

10  IF  (ITYP.EQ. 1)  Y2F=Y2F+ DELAF/ AO 


*  e  (de/dvbar) 

Y( 1 , 1) =E0 

«  a  bar ,  dabar/dvbar 
Y(2 , 1) =A0/A0 

*  t  (dt/dvbar) 

Y(3 , 1) =0.0 

*  nu  (dnu/dvbar) 

Y(4 , 1) =0 . 0 


*  Initial  Mean  Anomaly 

MAN0M0=DSQRT(MU/A0*»3) 

PI=DBLE(ACOS (- 1 . ) ) 
PI2=2 . 0*PI 
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o  o  o 


C 


C 

C 

C 

C 

C 

C 


C 


C 


C 


*  Initialize  ra*  &  rp* 

RABAR( 1) =Y(2 , 1) *  (1 . DO+EO) 
RABARI =RABAR ( 1 ) 
RABARF=Y2F* ( 1 . DO+EF) 
RPBAR( 1) =Y(2 , 1) * ( 1 . DO-EO) 
RPBARI =RPBAR ( 1 ) 
RPBARF=Y2F* ( 1 . DO-EF) 


»  PRINT  INITIAL  DATA  * 
a######*##**#####*#*## 


IF  ( IOUT 

.EQ.  1) 

THEN 

WRITE 

(9,*) 

'  Orbit  Data* 

WRITE 

(9,*) 

’  »»«**»****' 

WRITE 

(9 ,  *) 

WRITE 

(9,15) 

E0,A0,Y(2, 1) 

WRITE 

(9,16) 

EF.AF.Y2FCK 

WRITE 

(9,*) 

WRITE 

(9,*) 

WRITE 

(9,*) 

’  Primary  Body  Data’ 

WRITE 

(9.*) 

1  a***#**#*###****# ' 

WRITE 

(9,*) * 

WRITE 

(9,17) 

MU,  R 

WRITE 

(9.*) 

WRITE 

(9,*) 

WRITE 

(9,*) 

'  Transfer  Vehicle  Data’ 

WRITE 

(9,*) 

’  a*##*#####*##*#*#**## ' 

WRITE 

(9.*) 

WRITE 

(9,18) 

MO,  ISP 

WRITE 

(9,19) 

MDOT , MDOTBAR 

WRITE 

(9,*) 

WRITE 

(9,*) 

ENDIF 

IF  (ITYP 

.EQ. 1) 

THEN 

WRITE 

(*,») 

’  «»*»*»»»*«»»»*«»«***•’ 

WRITE 

(* , 20) 

Y2F 

WRITE 

(»,*) 

'  *»»»»»*»»*»»»*»»»*»*»’ 

WRIT" 

(9,*) 

’  *********************’ 

WRITE 

(9,20) 

Y2F 

WRITE 

(9,*) 

’  *#*#**##*#»*###*##**#  ’ 

ENDIF 

WRITE  (*,*) 
WRITE  (*,21) 
WRITE  (*,*) 
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CJ  o  u  o  o 


c 


IF  (IOUT.GE. 1)  THEN 
WRITE  (9,*) 

WRITE  (9,21) 

WRITE  (9,0 
ENDIF 


11  FORMAT 

12  FORMAT 

13  FORMAT 
♦ 

14  FORMAT 

15  FORMAT 
+ 

16  FORMAT 
+ 

17  FORMAT 
♦ 

18  FORMAT 
+ 

19  FORMAT 
+ 

20  FORMAT 

21  FORMAT 
+ 


(3X, ’Vehicle  Initial  Mass:  \F9.2,'  kg’) 

(3X, ’Propellant  Mass  Flow  Rate:  ’,E13.7,’  kg/sec’) 

(3X, ’Nondimensional  Propellant  Mass  Flow  Rate:  ’,E13.7, 
’  /sec’ ) 

(IX, ’system  Isp  currently  set  at  ’,F6.1,'  sec  ?') 

(3X, ’Initial  Orbit  Eccentricity:  ’ ,F6. 3, 6X, ’Initial  ’, 
’Orbit  Semimajor  Axis:  ’.F10.2,’  km  ( ’ ,E13 .7 , ’ ) ’ ) 

(3X, ’Final  Orbit  Eccentricity:  ’ ,F6.3,6X, ’Final  ’, 

’Orbit  Semimajor  Axis:  *,F10.2,'  km  ( ' ,E13 .7 , ’ ) ’ ) 

(3X, ’Gravitational  Parameter:  ’.E13.6,’  km  /sec  ’ ,6X, 
‘Radius:  ’ ,F9.3, ’  km’) 

(3X, ’Total  Initial  Mass:  \F9.2,’  kg’ ,6X, ’Specif ic  ’, 
’Impulse:  ’,F9.2,’  sec’) 

(3X, ’Propellant  Mass  Flow  Rate:  ’ ,E13 . 7 , ’kg/sec  (’, 

E13 . 7 , ’  /sec)’) 

(3X, ’AF/AO  =  ’ ,E13.7) 

(3X, ’Rev’ ,5X, ’Time  (sec) ’ ,4X, ’NO  (Rad) ' ,7X, ' V-bar ’ , 14X, 
’e’13X, ’a-bar’ 12X, ’ra-bar’ , 11X, ’rp-bar’) 


*  BEGIN  ITERATIONS  OF  STATE  VECTOR  * 
*«<hhhhhhh(***********#*«***#*»*#*»»* 


C 

C 

c 

c 

c 

c 

c 

c 


NREV= 1 
Y4NU=Y(4 , 1) 
VBAR=0 . DO 


*  Identify  Initial  Transfer  Procedure 

*  CASE  11  -  Increasing  a  and  e 

*  CASE  12  -  Increasing  a  and  Decreasing  e 

«  CASE  21  -  Decreasing  a  and  Increasing  e 

»  CASE  22  -  Decreasing  a  and  e 

FLAG=0 
FLAG1=- 1 

DIFFAsDABS (RABARF-RABARI ) 

DIFFP=DABS (RPBARF-RPBARI ) 

ICASE=21 

IF  ((RPBARF.GT.RPBARI) .OR. (DIFFP.LT. l.E-3))  THEN 
FLAG1=1 

IF  (DIFFP.LT. l.E-3)  FLAG=FLAG+ 1 
ICASE=11 

IF  (RABARF . LT . RABARI )  ICASE= ICASE+ 1 1 
IF  (DIFFA.LT. l.E-3)  GO  TO  70 
ENDIF 
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25  NXT=0 

icnt=o 

CALL  HAMING (NXT) 

IF  ((NXT.EQ.O) .AND. (H.LT.l.D-6))  THEN 
WRITE  (»,*)  *H  decreased  to  ’.H,’.’ 

WRITE  (*,*)  ’ Haming  still  refused  to  initialize.’ 

STOP 

ENDIF 

IF  (NXT.EQ.O)  THEN 
H=H/1.1 
IPNT=30 
GOTO  25 
ENDIF 

WRITE  (* , 130)  NREV, Y (3 , 1) , Y4NU , VBAR, Y (1,1) ,Y(2,1) , 

+  RABAR(l) .RPBAR(l) 

IF  (IOUT.GE. 1)  WRITE  (9,130)  NREV, Y(3 , 1) ,Y4NU, VBAR. Y( 1 , 1) , 

+  Y(2,l) ,RABAR(1) ,RPBAR(1) 

50  DO  100  1=1, I MAX 

CALL  HAMING (NXT) 

NREV= 1+INT (Y(4 ,NXT) /PI 2) 

Y4NU=Y (4 , NXT) -FLOAT (NREV- 1 ) *PI2 

IF  (ICASE.EQ.il)  DIFF=RABARF-RABAR(NXT) 

IF  (ICASE.EQ.22)  DIFF=RABAR(NXT) -RABARF 
IF  (ICASE.EQ. 12)  DIFF=RPBARF-RPBAR(NXT) 

IF  (ICASE.EQ. 21)  DIFF=RPBAR (NXT) -RPBARF 
IF  (DABS(DIFF) .LT. l.D-13)  GO  TO  110 
IF  (DIFF.GT.O.DO)  GO  TO  90 

WRITE  (*,») 

WRITE  (*,130)  VBAR, Y(l, NXT) ,Y(2, NXT) .RABAR(NXT) .RPBAR(NXT) 

IF  (IOOT.EQ. 1)  THEN 
WRITE  (9 , *) 

WRITE  (9,130)  VBAR, Y( 1 , NXT) ,Y(2,NXT) .RABAR(NXT) .RPBAR(NXT) 
ENDIF 

CALL  HAMING (NXT) 

WRITE  (»,*) 

WRITE  (*,130)  VBAR, Y ( 1 , NXT) ,Y(2, NXT) ,RABAR(NXT) .RPBAR(NXT) 
WRITE  (»,*) 

IF  (IOUT.EQ. 1)  THEN 
WRITE  (9 , ») 

WRITE  (9,130)  VBAR, Y ( l.NXT) ,Y(2, NXT) .RABAR(NXT) ,RPBAR( NXT) 
WRITE  (9,*) 
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END  IF 


NXT 1= NXT 
NXT=  1 

VBAR1=VBAR 
DO  55  J=  1 , 4 
Jl=5-J 

IF  (dCASE.EQ.il)  .OR.  (ICASE.EQ. 22))  X ( J 1 ) =RABAR(NXT1) 
IF  ((ICASE.EQ. 12) .OR. (ICASE.EQ. 21))  X(J1) =RPBAR(NXT1) 
VAR1 ( J 1 ) =VBAR1 
VAR2 (J1 ) =Y( 1 ,NXT1) 

VAR3 ( J 1 )  =  Y  ( 2  , NXT 1 ) 

VAR4 ( J1 ) =Y(3 , NXT1) 

VAR5 ( J 1 ) = Y (4 ,NXT1 ) 

VBAR1=VBAR1-H 
NXT1=NXT1- 1 
IF  (NXT1.EQ.0)  NXT 1=4 
CONTINUE 

IF  ((ICASE.EQ.il) .OR. (ICASE.EQ. 22))  XBAR=RABARF 
IF  ((ICASE.EQ. 12) .OR. (ICASE.EQ. 21))  XBAR=RPBARF 
CALL  INTERP3  (XBAR.X, VAR1 ,P) 

VBAR=P 

CALL  INTERP3  (XBAR.X, VAR2 ,P> 

Y ( 1 ,NXT) =P 

CALL  INTERP3  (XBAR.X, VAR3 ,P) 

Y(2 ,NXT) =P 

CALL  INTERP3  (XBAR.X, VAR4 ,P) 

Y(3,NXT)=P 

CALL  INTERP3  (XBAR.X, VAR5 ,P) 

Y (4 ,NXT) =P 

IF  (dCASE.EQ.il)  .OR.  (ICASE.EQ.  22))  RABAR(NXT)  =RABARF 
IF  ((ICASE.EQ. 12) .OR. (ICASE.EQ. 21))  RPBAR(NXT) =RPBARF 
NREV= 1+INT (Y(4 ,NXT) /PI 2) 

Y4NU=Y (4 ,NXT) -FLOAT (NREV- 1 } *PI2 

FLAG=FLAG+ 1 

IF  (FLAG.EQ.2)  GOTO  80 
ICASE=12 

IF  (FLAG1.LT. 0)  THEN 
ICASE= 1 1 

IF  (RABARF.LT. RABARI)  ICASE=ICASE+11 
IF  (DIFFA.LT. l.E-3)  GO  TO  70 
END  IF 

GO  TO  25 


WRITE  (*,130)  NREV , Y(3 ,NXT) , Y4NU , VBAR, Y ( 1 ,NXT) , Y (2 ,NXT) , 
RABAR(NXT) , RPBAR(NXT) 

WRITE  (*,*) 

IF  (IOUT.GE. 1)  THEN 


o  o  o  o  o  o  o  o  a  non  non  n  n  n  n  n  n 


WRITE  (9,130)  NREV, Y(3 ,  NXT) , Y4NU, VBAR, Y( 1 ,NXT) , Y (2 ,NXT) , 
+  RABAR(NXT) .RPBAR(NXT) 

WRITE  (9 ,  *) 

END  IF 
GO  TO  150 


90  IF  (IPNT.EQ.O)  GO  TO  100 

IF  ( (I/IPNT»IPNT) .EQ. I)  THEN 

WRITE  (#,130)  NREV,Y(3,NXT) ,Y4NU,VBAR,Y(1 ,NXT) ,Y(2,NXT) , 

♦  RABAR(NXT) , RPBAR(NXT) 

IF  (IODT.EQ. 1)  WRITE  (9,130)  NREV,Y(3,NXT) .Y4NU.VBAR, 

+  Y(1,NXT) ,Y(2,NXT) , RABAR(NXT) .RPBAR(NXT) 

ENDIF 


100  CONTINUE 


WRITE  (»,#)  'Maximum  number  of  iterations  reached.' 
WRITE  (»,*)  'Program  Terminated’ 

GOTO  170 


110  Y( 1 , 1) SY( 1 ,NXT) 

Y(2 , 1) =Y(2 ,NXT) 

Y(3 , 1 ) =Y(3 ,NXT) 

Y(4 , 1) =Y(4 ,NXT) 
RABAR(l) =RABAR(NXT) 
RPBAR( 1) =RPBAR(NXT) 

GO  TO  70 


130  FORMAT  (3X, 13 ,3X,E14 . 7 ,3X, F6 . 3 , 3X , E14 . 7 , 3X.E14 . 7 ,3X, 
♦  E14.7,3X,E14.7,3X,E14.7) 

135  FORMAT  (3X. 'de/dVbar  =  ’.E20.13) 

150  CONTINUE 


»  TRANSFER  TOTAL  DELTA-V  AND  TIME  CALCULATION  » 

I****#*##*#*###**###*##*##*#*****##**#**##*#*### 

DELV=VBAR»AO«MANOMO 
DEL VM=DELV* 1000 . 

T0F= 1 . /MDOTBAR# ( 1 . -DEXP ( -DELV/ (9.8 1D-3»ISP) ) ) /3600 . /24 . 
TOFY=TOF/365 . 25 


#  Comparison  Data  (Circular  Orbits  Only)  » 
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o  o  o  o  o  o  o  o 


c 

C 


###»***»*##******#*##* 


IF  ((E0.LT.1.E-6) .AMD. (EF.LT.l.E-6))  GOTO  170 


*  Spiral  Transfer 

VCS 1 =DSQRT ( MO/ AO ) 

VCS2sDSQRT (MU/A0/Y2F) 

DELVSPaVCSl-VCS2 
DEL VSMaDELVSP* 1000. 

VSPBARaDELVSP/AO/MANOMO 

T0FSP= 1 . /MDOTBAR* ( 1 . -DEXP (-DELVSP/ (9 . 81D-3*ISP) ) ) /3600 . /24 . 
TOFSPYaTOFSP/365 . 25 


*  Hohmann  Transfer 


C 

C 


C 


V1=DSQRT(2»MU* ( 1 . /AO- 1 . / <A0+A0»Y2F) ) ) 
V2=DSQRT (2*MU* ( 1 . /A0/Y2F- 1 . / (A0+A0»Y2F) ) ) 
DELVHM= (V1-VCS1 ) + (VCS2-V2) 

DEL VHMMaDELVHM» 1000 . 

VHMB AR= DELVHM/ AO / MANOMO 


WRITE 

(*,») 

WRITE 

(» , 160) 

DELV.DELVM 

WRITE 

(« , 165) 

TOF.TOFY 

WRITE 

(«,«) 

WRITE 

(«  ,«) 

WRITE 

(*,*)  ’ 

Comparison  Transfer  Data’ 

WRITE 

(*,»)  ’ 

************************ ’ 

WRITE 

(»,*) 

WRITE 

(*,*)  ’ 

Spiral  Transfer' 

WRITE 

(*,«) 

WRITE 

(* ,  160) 

DELVSP .DELVSM 

WRITE 

(» , 165) 

TOFSP.TOFSPY 

WRITE 

(*,») 

WRITE 

(«,«)  ’ 

Hohmann  Transfer’ 

WRITE 

(*,») 

WRITE 

(*,160) 

DELVHM .DELVHMM 

WRITE 

(*  ,*) 

WRITE 

(*  ,*) 

WRITE 

(*  ,*) 

IF  (IODT.GE. 1)  THEN 
WRITE  (9 , #) 

WRITE  (9,160)  DELV.DELVM 
WRITE  (9,165)  TOF.TOFY 
WRITE  (9 ,  *) 

WRITE  (9,») 

WRITE  (9,*)  ’  Comparison  Transfer  Data’ 
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no  on  on 


WRITE 

(9,*)  ’ 

#»»*#»###########****###’ 

WRITE 

(9,*) 

WRITE 

(9 , »)  ’ 

Spiral  Transfer’ 

WRITE 

(9 ,  *) 

WRITE 

(9,160) 

DELVSP.DELVSM 

WRITE 

(9,165) 

TOFSP.TOFSPY 

WRITE 

(9,*) 

WRITE 

(9 , »)  ’ 

Hohmann  Transfer’ 

WRITE 

(9 , « ) 

WRITE 

(9,160) 

DELVHM , DEL VHMM 

WRITE 

(9,*) 

WRITE 

(9,*) 

WRITE 

(9 , ») 

ENDIF 

’  (ITYP 

.EQ.  1) 

THEN 

WRITE 

( 10 , *) 

Y2F.VBAR 

WRITE 

(11,*) 

Y2F , VSPBAR 

WRITE 

(12, *) 

Y2F.VHMBAR 

END  IF 


C 

•  160  FORMAT  (3X,’Total  Transfer  Delta-V:  ’ ,F9. 

+  F8.2,  ’  m/s) ’) 

165  FORMAT  (3X,’Total  Transfer  Time:  ’,F9.3,’ 
+  F6.2,  ’  yr) *) 

C 

C 

•  DIFFYF=Y2FCK-Y2F 

IF  (DABS(DIFFYF) .GT.l.E-13)  GOTO  10 


170  CONTINUE 


IF  (I0UT.GE.1)  CLOSE  (9) 
IF  (ITYP.EQ. 1)  THEN 
CLOSE  (10) 

CLOSE  (11) 

CLOSE  (12) 

ENDIF 


C 

C 

C 

C 

C 

C 

C 


STOP 

END 


##*«#«*##»»**#*#*#*## 
»  SUBROUTINE  RHS 

»**»»*»•»**«»»»«****» 

SUBROUTINE  RHS(NXT) 


4 , *  km/s  ( ’ , 
days  ( ’ , 


*»»»*****«*« 

* 

**«*»*»***«* 


C-l.  11 


o  o  o  non  can  <n  ana  a  o  o  o 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  MANOMO , MDOTBAR 

COMMON  /HAM/  VBAR, Y(42 ,4) , F (42 , 4) , ERR(42) , N, H.MODE 
•  COMMON  /RH/  RABAR(4) .RPBAR (4) .ICASE, I CNT, MANOMO, AO .MDOTBAR, 

♦  ACCEL0.PI2 


*  Specify  ra  =  constant  (ID=0)  or  rp  =  constant  (ID=1) 

IF  ((ICASE.EQ.il) .OR. (ICASE. EQ. 22))  ID=1 
IF  ((ICASE. EQ. 12). OR. (ICASE. EQ. 21))  ID=0 

CALL  INTERPO (Y ( 1 , NXT) ,DELA , DELE , ID , ICNT) 

ICNT=ICNT+ 1 


*  Define  Derivatives 

IF  ((ICASE. EQ. 11) .OR. (ICASE. EQ. 21))  SIGN=1.0 
IF  ( (ICASE. EQ. 12) .OR. (ICASE. EQ. 22))  SIGN=-1.0 

F ( 1 ,NXT) =Y (2 ,NXT) ** . 5*SIGN*DELE/PI2 

F (2 , NXT) =Y(2 , NXT) **1 . 5*SIGN*DELA/PI2 

F (3 ,NXT) = ( 1 . -MDOTBAR*Y(3 ,NXT) ) * MANOMO* AO /ACCELO 

Y4NU1=Y (4 ,NXT) -FLOAT (INT (Y(4 ,NXT) /PI2) ) 

F (4 ,NXT) =MANOMO/ (Y(2 ,NXT) * (1 . -Y( 1 ,NXT) **2) ) ** 1 . 5* 
+  ( 1 . +Y( 1 ,NXT) *DCOS (Y4NU1) ) **2*F (3 ,NXT) 

*  Define  Changes  in  ra-bar  and  rp-bar 

RABAR (NXT) =Y ( 2 , NXT) * (1 . +Y(1 , NXT) ) 

RPBAR (NXT) =Y(2 ,NXT) * ( 1 . -Y( 1 , NXT) ) 

RETURN 

END 


SUBROUTINE  HAMING(NXT) 

C  VERSION  OF  11/20/1987 
C  PURPOSE 

C  HAMING  IS  AN  ORDINARY  DIFFERENTIAL  EQUATIONS  INTEGRATOR 
C  IT  IS  A  FOURTH  ORDER  PREDICTOR-CORRECTOR  ALGORITHM  WHICH 
C  MEANS  THAT  IT  CARRIES  THE  LAST  FOUR  VALUES  OF  THE  STATE 

C  VECTOR,  AND  EXTRAPOLATES  THESE  VALUES  TO  OBTAIN  A  PREDICTED 

C  NEXT  VALUE  (THE  PREDICTION  STEP)  AND  EVALUATES  THE  EQUATIONS 
C  OF  MOTION  AT  THE  PREDICTED  POINT,  AND  THEN  CORRECTS  THE 

C  EXTRAPOLATED  POINT  USING  A  HIGHER  ORDER  POLYNOMIAL  (THE 

C  CORRECTION  STEP) . 

C  INPUT 

C  NXT  --  IN  THE  CALL  SPECIFIES  WHICH  OF  THE  FOUR  VALUES  OF 
C  THE  STATE  VECTOR  IS  THE  CURRENT  ONE.  NXT  IS  UPDATED 

C  BY  HAMING  AUTOMATICALLY,  BUT  MUST  BE  SET  TO  ZERO  ON 


C-l. 12 


o  o  o  o 


SET  TO  REASONABLE  VALUE 


THE  FIRST  CALL. 

CALL  ROUTINES 
RHS(NXT) 

REFERENCES 

WILLIAM  WEISEL 
PROGRAMMER 

RODNEY  D.  BAIN 
PROGRAM  MODIFICATIONS 
NONE 

COMMENTS 

TOL  —  IS  HAMING’S  START  UP  TOLERANCE 

AS  NECESSARY 
THE  COMMON  BLOCK  CONTAINS: 

VBAR  --  IS  THE  INDEPENDENT  VARIABLE  (OFTEN  TIME) 

Y (42 ,4)  --  IS  THE  STATE  VECTOR,  4  COPIES  OF  IT,  WITH  NXT  POINTING 
POINTING  TO  THE  CURRENT  ONE,  THE  LIMIT  OF  42  EQUATIONS 
OF  MOTION  CAN  BE  CHANGED. 

ARE  THE  EQUATIONS  OF  MOTION  EVALUATED  AT  THE  SAKE  TIMES 
AS  THE  STATE  VECTOR  Y  ...  IT  IS  THE  JOB  OF  SUBROUTINE 
RHS  TO  CALCULATE  THESE. 

IS  AN  ESTIMATE  OF  THE  ONE-STEP  INTEGRATION  ERROR 
IS  THE  NUMBER  OF  ODES  ...  LIMIT  IS  42  UNLESS  YOU  CHANGE 
THE  COMMON  BLOCK 

IS  THE  TIMESTEP  ...  ONE  CALL  TO  HAMING  INCREMENTS  X  BY  H 
IS  ZERO  FOR  EOM  ONLY,  1 
THE  USER  MJST  SUPPLY  A  MAIN  PROGRAM, 

EVALUATES  THE  EQUATIONS  OF  MOTION. 


F  (42 , 4)  - 


ERR (42) 


H 

MODE 


FOR  EOM  AND  EOV 

AND  THE  SUBROUTINE  RHS (NXT) 


WHICH 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  /HAM/  VBAR,Y(42,4) ,F(42,4) ,ERR(42) ,N,H,MODE 
DATA  ZERO , ONE , TWO , THREE , FOUR/O . DO , 1 . DO , 2 . DO , 3 . DO . 4 . DO/ 
TOL= 1 . D- 12 


CHECK  IF  THIS  IS  THE  FIRST  CALL  . . .  HAMING  (LIKE  ALL  PREDICTOR- 
CORRECTORS)  NEEDS  'PREVIOUS'  VALUES 


IF (NXT)  190,10,200 
C 

C  IT  IS  A  PICARD  INTERATION  (SLOW  AND  EXPENSIVE)  TO  STEP  BACKWARDS 

C  IN  TIME  THREE  STEPS  TO  GET  THE  4  PREVIOUS  POINTS.  A  SUCCESSFUL 

C  STARTUP  RETURNS  NXT=1,  AND  TIME  HAS  NOT  BEEN  INCREMENTED.  IF 

C  STARTUP  FAILS,  NXT  WILL  BE  RETURNED  AS  ZERO. 

C 

10  XO=VBAR 
HH=H/TWO 
CALL  RHS ( 1 ) 

DO  40  L=2 , 4 
VBAR=VBAR+HH 
DO  20  1  =  1, N 

Y(I,L)=Y(I,L-1) +  HH*F (I , L- 1 ) 

20  CONTINUE 

CALL  RHS (L) 
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o  o  o  o 


• 

30 

VBAR=VBAR+HH 

DO  30  1=1, N 

Y(I,L)=Y(I,L-1)+H*F(I,L) 

CONTINUE 

40 

CALL  RHS(L) 

CONTINUE 

50 

JSW=-10 

ISW=1 

• 

DO  120  1=1, N 

HH=Y(I , 1) +H* (9 . DO»F(I , 1) ♦ 19 .DO»F (I , 2) -5 .DO*F (I ,3) 

70 

1  +F (I , 4) ) /24 . DO 

IF(DABS(HH-Y(I ,2) ) .LT.TOL)  GOTO  70 

ISW  =0 

Y(I ,2) =HH 

• 

HH=Y(I , 1) +H* (F(I , 1 ) +FOUR«F(I , 2) +F(I ,3) ) /THREE 

IF (DABS (HH-Y(I ,3)) .LT.TOL)  GOTO  90 

90 

ISW=0 

Y(I ,3) =HH 

• 

HH= Y (1,1) +H* (THREE»F ( I , 1 ) +9 . DO»F ( 1 , 2) +9 . DO*F ( 1 , 3) 

1  +THREE*F ( I , 4) ) /8 . DO 

IF (DABS (HH-Y(I ,4)) .LT.TOL)  GOTO  110 

ISW=0 

110 

Y(I ,4) =HH 

120 

CONTINUE 

• 

VBAR=XO 

DO  130  L=2 , 4 

VBAR=VBAR+H 

CALL  RHS(L) 

130 

CONTINUE 

140 

IF(ISW)  140,140.150 

JSW=JSW+1 

150 

IF(JSW)  50,280,280 

VBAR=XO 

• 

ISW=1 

160 

JSW=1 

DO  160  1=1, N 

ERR(I) =ZERO 

CONTINUE 

• 

NXT=  1 

GOTO  280 

C 

C  A  CALL  TO  HAMING  WITH  NXT=-NXT,  AFTER  A  SUCCESSFUL  STARTUP, 

C  WILL  TURN  OFF  THE  SECOND  EVALUATION  OF  THE  EQUATIONS  OF  MOTION 
C  FOLLOWING  THE  CORRECTOR  STEP.  IN  SYSTEMS  WHERE  THE  EQUATIONS  OF 
C  MOTION  ARE  VERY  EXPENSIVE,  THIS  CAN  HALVE  YOUR  RUN  TIME. 

C 

190  JSW=2 

NXT=IABS (NXT) 

THIS  IS  THE  PREDICTOR-CORRECTOR  ALGORITHM  . . .  FIRST  THE  INDICES 
ARE  PERMUTED 


C-l.  14 


o  o  o  o  o  o  o  o  o  o  o  o  oooo 


200  VBAR=VBAR+H 

NP1=M0D (NXT ,  4)  ♦  1 
GOTO  (210,230) ,ISW 

•  210  GOTO  (270,270,270,220) , NXT 

220  ISW=2 

230  NM2=M0D(NP1 , 4) ♦ 1 
NM1-M0D (NM2 ,  4)  *  1 
NP0=M0D (NM1 , 4) *  1 

. . .  THEN  THE  PREDICTOR  PART  IS  RUN  TO  FIND  AN  EXTRAPOLATED  VALUE 
OF  THE  STATE  VECTOR  AT  THE  NEW  TIME  ... 

DO  240  1=1, N 

F ( I , NM2) =Y(I ,NP1 ) +FOUR*H* (TWO#F ( I ,NP0) -F (I ,NM1 ) 

1  +TW0#F ( I , NM2 ) ) /THREE 

Y(I ,NP1) =F(I , NM2) -0 . 9256 19835D0»ERR ( I ) 

240  CONTINUE 

THE  EQUATIONS  OF  MOTION  ARE  EVALUATED  AT  THE  EXTRAPOLATED  VALUE 
OF  THE  STATE  VECTOR  . . . 

CALL  RHS(NPl) 

AND  THE  CORRECTOR  ALGORITHM  IS  USED  TO  ADD  THIS  NEW  INFORMATION 
AND  OBTAIN  A  BETTER  VALUE  OF  THE  NEW  STATE  VECTOR  . . . 

DO  250  1=1, N 

Y(I,NP1)=(9.D0«Y(I ,NPO) -Y(I , NM2) +THREE*H* (F (I ,NP1 ) 

1  +TWO»F (I ,NPO) -F ( I ,NM1) ) ) /8 . DO 

ERR ( I ) =F ( I , NM2 ) - Y ( I , NP 1 ) 

Y ( I ,NP1) =Y(I ,NP1)+0.0743801653D0»ERR(I) 

250  CONTINUE 

GOTO  (260,270) ,JSW 

FINALLY,  THE  EQUATIONS  OF  MOTION  ARE  RE-EVALUATED  AT  THE  BETTER 
VALUE  OF  THE  STATE  VECTOR  ...  THIS  CAN  BE  SUPPRESSED. 


260 

270 

280 


CALL  RHS(NPl) 

NXT=NP1 

RETURN 

END 


»»«»««»*«»««**»»««**»***»«»«*«*»* 


«  SUBROUTINE  INTERP3  * 

*  * 

*  The  following  subroutine  interpolates  between  four  data  » 

*  points  using  a  third  order  polynomial.  * 


»»»»»»«*••»»«*»«»»»*«*»*»»*»*»*»* 
SUBROUTINE  INTERP3  (XBAR.X.D.P) 
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IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  X ( 4 ) ,P1 (3) ,D(4) 

C 

DO  300  K= 1 , 3 

DO  300  J= (K+ 1)  ,4 
JI=4-J+ (K+l) 

D(J1)=(D(J1)-D(J1-1))/(X(J1) -X(Jl-K) ) 

300  CONTINUE 
C 

Z=XBAR-X( I) 

PI ( 1) SD ( 1) +D (2) *Z 
DO  310  J=2 , 3 
Z=Z# (XBAR-X(J) ) 

P1(J)=P1(J-1)+D(J+1)»Z 
310  CONTINUE 
C 

P=P1 (3) 

C 

RETURN 

END 

C 

C 

C 

C 

C  **«**«*»#**##»««*«»»**»»»»»»»** 

C  «  SUBROUTINE  INTERPO  » 

C  ###***»#»*#*»#***#»#»#»*#*«**## 

C 

c  *»»#*##*######*###*««*»***»*»#* 

C  *  The  following  subroutine  utilizes  the  Newton  formula  * 

C  *  to  interpolate/extrapolate  from  the  given  data  sets  » 

C  *  utilizing  a  4th  degree  interpolating  polynomial.  » 

C  a#*#*###########**######**##### 

c 

c 

SUBROUTINE  INTERPO (EBAR , DA , DE , ID , ICNT) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E ( 100) .LAMDA(IOO) .DELA(IOO) .DELEUOO) 

C 

C 

C  »  OPEN  DATA  FILES  * 

C  »*###*###***«****## 

c 

C  »  Constant  ra  Data 

C 

IF  ((ID.EQ.O) .AND. (ICNT.EQ.O))  THEN 

OPEN  (UNIT= 15 ,  STATUS  =  ’OLD’ .FILE  =  ’ ALAM.DAT’ ) 

OPEN  (UNIT= 18 ,  STATUS  =  ’OLD’, FILE  =  ’ADELA.DAT’ ) 

OPEN  (UNIT® 17 ,  STATUS  =  ’OLD’, FILE  =  ’ADELE.DAT’ ) 

IDATA=99 

ENDIF 
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*  Constant  rp  Data 

IF  ( (ID.EQ. 1) .AND. (ICNT.EQ.O) )  THEN 
OPEN  (UNIT= 15 ,  STATUS  =  ’OLD’, FILE  =  ’PLAM.DAT’ ) 

OPEN  (UNIT= 16 ,  STATUS  =  'OLD '.FILE  =  'PDELA.DAT' ) 

OPEN  (UNIT=17,  STATUS  =  'OLD' .FILE  =  'PDELE.DAT') 

IDATA=98 
END  IF 


C 


C 

C 

C 

C 

C 


C 


C 


C 


C 

C 

c 

c 

c 

c 

c 


c 


IF  (ICNT.EQ.O)  THEN 

BEWIND  (15) 

REWIND  (16) 

REWIND  (17) 


*  INPUT  INITIAL  DATA  « 
a#*###*##**##********* 

DO  5.  1=1 , IDATA 

READ  (15.*)  E(I) ,LAMDA(I) 
READ  (16.*)  E(I).DELA(I) 
READ  (17.*)  E(I) .DELE (I) 
CONTINUE 

IDEG=4 

CLOSE  (15) 

CLOSE  (16) 

CLOSE  (17) 

END  IF 


*  EVALUATE  FUNCTION  AT  SPECIFIED  VALUES  OF  ECCENTRICITY  * 
********************************************************* 


*  Check  location  of  data  point  * 

DO  30  1=1, IDATA 
DIFF=EBAR-E(I) 

IF  (ABS(DIFF) .LE. l.E-15)  THEN 
LAM=LAMDA(I) 

DA=DELA(I) 

DE=DELE(I) 

GO  TO  40 
END  IF 

IF  ( (I.EQ. 1) .AND. (DIFF.LT.O. ) )  THEN 
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C 


C 


C 


C 

C 


IMIN=  1 

IMAX=IMIN+IDEG 
GO  TO  35 
END  IF 

IF  ((I.EQ.IDATA) .AND. (DIFF.GT.O.))  THEN 
I  MAXIDATA 
I MI N= I MAX- I DEG 
GO  TO  35 
ENDIF 

I LEFT = I 

IF  (DIFF.LT.O. )  GO  TO  31 
30  CONTINUE 


31  IDEGP1=IDEG+ 1 
IDEGL=IDEGPl/2 
I MI N= I LEFT  +1-1 DEGL 
IMAX=IMIN+IDEG 

IF  (IMIN.LT. 1)  THEN 
IMIN=1 

IMAX=IMIN+IDEG 

ENDIF 

IF  (IMAX.GT.IDATA)  THEN 
IMAX=IDATA 
I MI N= I MAX- I DEG 
ENDIF 


*  Evaluate  Polynomial  * 


C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


35  CALL  INTEBP(IDEG.IMIN,IMAX,EBAR,E,LAMDA, DELA.DELE.LAM.DA.DE) 

40  CONTINUE 
RETURN 
END 


«*«**»*»*»»»»*»»»*»**»*»»«»*»*** 
»  SUBROUTINE  INTERP  » 

#»*»#*##»*»########*##»*#######* 

*»###»##*#*#*##»#*###*#**#»**#*# 

*  The  following  algorithm  evaluates  an  interpolating  * 

*  nth  degree  polynomial  * 

»»»»»»»»»*«»«*»»»*«»»»«»»»»»»«»» 
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c 

SUBROUTINE  INTERP ( IDEG . I MIN , IMAX , EBAR , E , LAMDA , DELA . DELE , LAM, DA , BE) 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E(100) , LAMDA ( 100) .DELA(IOO) .DELE(IOO) 

DIMENSION  X(100) ,D(100) ,P(IOO) 

C 

IDEGP1=IDEG+ 1 
C 

DO  200  1*1,3 
C 

DO  150  J-IMIN.IMAX 
K=J-IMIN+ 1 

IF  (I.EQ.l)  D (K) =LAMDA(J) 

IF  (I.EQ.2)  D (K) =DELA(J) 

IF  (I.EQ.3)  D (K) =DELE (J) 

X(K)=E(J) 

150  CONTINUE 
C 

DO  155  K= 1 , IDEG 

DO  155  J= (K+l) , IDEGP1 
J1=IDEGP1-J+ (X+l) 

D(J1)S(D(J1)-D(J1-1))/(X(J1) -X(Jl-K) ) 

155  CONTINUE 

C 

Z=EBAR-X( 1) 

P(l) =D(1) +D(2) *Z 
C 

DO  160  J=2 , IDEG 
Z=Z* (EBAR-X(J) ) 

P(J)=P(J-1)+D(J+1)»Z 
160  CONTINUE 
C 

IF  (I.EQ.l)  LAM=P (IDEG) 

IF  (I.EQ.2)  DA=P(IDEG) 

IF  (I.EQ.3)  DE=P(IDEG) 

C 

200  CONTINUE 
C 

RETURN 

END 
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Appendix  C-2: 


Program  TRAHSALT 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


**#########*»#*#########**###*#*# 


«  * 

*  Program  TBANSALT  * 

»  « 

*  Captain  Gregory  Beeker  « 

*  Air  Force  Institute  of  Technology  » 

*  July  1988  * 

«  * 


#*########»*##*#*»#*#*###*#*##*## 

The  following  program  was  written  to  solve  the  long  timescale 
problem  to  determine  the  total  transfer  time  of  a  spacecraft 
traveling  between  two  orbits.  The  program  takes  a  transfer 
vehicle  from  an  initial  orbit  to  a  final  orbit  though  many 
revolutions.  A  combination  of  the  fast  timescale  changes  in  the 
orbital  elements  for  constant  distance  to  perigee  (rp)  and  apogee 
(ra)  over  two  revolutions  is  implemented.  Subroutine  Haming 
is  incorporated  to  solve  the  differential  equations  for  the 
nondimensional  changes  in  eccentricity  and  semimajor  axis 
(dabar/dVbar  and  de/dVbar) .  In  addition,  Program  INTERPO,  which 
provided  the  final  solutions  to  the  fast  timescale  problem,  is 
also  incorporated  as  a  subroutine  to  provide  the  values  of  delta 
a»  and  delta  e*  as  functions  of  eccentricity. 


IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  I SP , MU , MANOMO , MANOM , MDOTBAR , MDOT , MO 
DIMENSION  X(4) ,VAR1(4) ,VAR2(4) ,VAR3(4) ,VAR4(4) 

COMMON  /HAM/  VBAR, Y(42 , 4) , F (42 , 4) ,ERR(42) ,N,H , MODE 
COMMON  /RH/  RABAR(4) ,RPBAR(4> , ICASE, ICNT.DIFF 


C 

C 

C 

C 

C 

C 

C 

C 


»  INITIALIZATION  OF  PARAMETERS  FOR  HAMING  » 
••a#***###*#**###*#**##*****#**###*###**#*# 

N=2 

M0DE=0 


«  INITIAL  DATA  ENTRY  * 

»»»«»*«»**«»********«* 

#  Primary  Body  Parameters 

WRITE  («,»)  ’Do  you  wish  to  change  the  value  of  the’ 
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+ 


WHITE  (»,»)  ’gravitational  parameter  and  radius  of  the’ 

WRITE  (*,*)  'primary  body  (currently  set  to  those  of’ 

WRITE  (*,*)  'the  Earth)?  (0  -  yes;  1  -  no)’ 

+  READ  (*,»)  IP 

C 

IF  (IP.EQ.O)  THEN 

WRITE  (*,*)  ’Enter  the  new  values  of  mu  (km3/sec2)  and  r  (km).’ 
READ  (*.*>  MU,R 
END  IF 

#  IF  (IP.EQ.l)  THEN 

MO-3 . 986012D5 
R=8378. 145D0 
END  IF 

*  Orbit  Data 

WRITE  (*,#)  ’Specify  the  type  of  program  run  (0  or  1).’ 

WRITE  (*,*) 

WRITE  (*,*)  ’  0  -  Single  Transfer’ 

WRITE  (»,»)  ’  1  -  Multiple  Transfers  (Includes  plot’ 

WRITE  (»,«)  ’  data  file  of  af/aO  vs  delta-V  total)’ 

READ  (»,*)  ITYP 

WRITE  (»,«)  'Enter  the  eccentricity  and  3emimajor  axis  (km)’ 

WRITE  (*,*)  'of  the  initial  orbit.’ 

READ  (»,»)  E0,A0 

2  IF  (AO.LE.R)  THEN 

WRITE(»,»)  'Initial  orbit  intersects  the  surface  of  the’ 

WRITE (*,*)  'primary  body.  Please  select  a  new  value.’ 

READ  (*,«)  AO 
GOTO  2 
ENDIF 

WRITE  (»,*)  ’Enter  the  eccentricity  and  semimajor  axis  (km)’ 

WRITE  (»,*)  'of  the  final  orbit.' 

READ  (*,*)  EF , AF 

3  IF  (AF.LE.R)  THEN 

WRITE(»,»)  ’Final  orbit  intersects  the  surface  of  the’ 
WRITE(*,»)  ’primary  body.  Please  select  a  new  value.' 

READ  (*,»)  AF 
GOTO  3 
ENDIF 

IF  (ITYP.EQ. 1)  THEN 

WRITE  (*,*)  ’Specify  the  number  of  transfers  between’ 

WRITE  (»,*)  ’the  initial  and  final  orbits.’ 

READ  (*,»)  NTRANS 
ENDIF 

«  Vehicle  Parameters 
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C 


M0=2000 
MD0TS4 . OD-4 

ft  MDOTBAR=MDOT/MO 

ISP=5000 

C 

WRITE  (*,*)  ’Do  you  wish  to  change  the  value  of  the  vehicle  mass’ 
WRITE  (*,»)  'and  mass  flow  rate  parameters?  (0  -  yes;  1  -  no)’ 
WRITE  (*,*) 

•  WRITE  (*,*)  ’Current  Settings:’ 

WRITE  (*.11)  MO 

WRITE  (* , 12)  MDOT 
WRITE  (*,13)  MDOTBAR 
READ  (*,*)  IM 
C 

•  IF  (IM.EQ.O)  THEN 

WRITE  (*,*)  'Enter  the  new  values  of  the  mass  flow  rate  (kg/s)’ 
WRITE  (*,»)  ‘and  initial  vehicle  mass  (kg).’ 

READ  (*.»)  MDOT , MO 
MDOTBAR=MDOT/MO 
WRITE  (*.13)  MDOTBAR 

•  ENDIF 
C 

WRITE  (*,*)  ’Do  you  wish  to  change  the  value  of  propulsion’ 

WRITE  (*.14)  ISP 

WRITE  (*.«)  ’ (0  -  yes;  1  -  no) ’ 

READ (• . *) IISP 

•  C 

IF  (IISP.EQ.O)  THEN 

WRITE  (*,*)  ’Enter  the  new  value  of  Isp.’ 

READ  (*,*)  ISP 
ENDIF 


C 


*  Integration  Data 


WRITE  (*,*)  'Enter  the 
WRITE  (*,*)  'variable, 
READ  (»,*)  H 
WRITE  (*,*)  'Enter  the 
WRITE  (*,*)  'performed 
READ  (»,«)  IMAX 


step  size  of  the  independent’ 

V-bar. ’ 

maximum  number  of  iterations  to  be' 
by  Haming  during  each  transfer. ’ 


*  Output  Parameters 

IF  (ITYP.EQ. 1)  THEN 
I0(JT=2 

OPEN  (UNIT  =  10.  STATUS  =  ’NEW’,  FILE  =  ’EDVALT’) 

OPEN  (UNIT  =  11,  STATUS  =  ’NEW’,  FILE  =  ’ SDVALT’ ) 

OPEN  (UNIT  =  12,  STATUS  =  ’NEW’,  FILE  =  ’HDVALT’ ) 

GO  TO  8 
ENDIF 
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WRITE  (*,*)  'Enter  output  code.’ 

WRITE  («.*)  '(0  for  screen,  1  for  file  and  screen)’ 

READ  (*,»)  I OUT 
C 

8  IF  (IOUT.GE. 1)  THEN 

WRITE  («,*)  ’Specify  the  name  of  the  data  output  file. 
READ  (* , ’ (A) ’ )  TRANSOUT 

OPEN  (UNIT  =  9,  STATUS  =  ’NEW’,  FILE  =  TRANSOUT) 

ENDIF 

C 

IPNT=0 

IF  (IOUT.LE. 1)  THEN 

WRITE  (*,*)  'Enter  the  number  of  steps  between  each’ 
WRITE  (*,»)  'data  printout.’ 

READ  (»,»)  IPNT 
ENDIF 
C 
C 
C 

C  »  INITIALIZE  ORBIT  PARAMETERS  (NONDIMENSIONAL)  * 

C  ***#*###*####****»**#*##**##*»##*#**###***##**#»# 

c 

C  *  State  Vector 

C 

Y2F=AF/A0 

Y2FCK=Y2F 

C 

IF  (ITYP.EQ. I)  THEN 
DELAF= (AF-AO) /NTRANS 
Y2F= 1 . DO 
ENDIF 
C 

10  IF  (ITYP.EQ. 1)  Y2F=Y2F+DELAF/A0 
C 

Y( 1 , 1) =E0 
Y ( 2 , 1 ) =A0/A0 


»  Initial  Mean  Anomaly 
MANOMO=DSQRT(MU/AO**3) 


#  Initialize  ra*  &  rp* 

RABAR ( 1 ) =Y(2 , 1 ) » ( I .DO+EO) 
RABARI =RABAR ( 1 ) 
RABARF=Y2F» ( 1 . DO+EF) 
RPBAR( 1)=Y(2,1)»(1. DO-EO) 
RPBARI=RPBAR( 1) 
RPBARF=Y2F» ( 1 . DO-EF) 
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«  Identify  Orbit  Type  « 
*********************** 

*  CASE  11  -  Increasing  a  and  e 

*  CASE  12  -  Increasing  a  and  Decreasing  e 

*  CASE  21  -  Decreasing  a  and  Increasing  e 

»  CASE  22  -  Decreasing  a  and  e 

ICASE=1 1 

IF  (AO.GT.AF)  ICASE=ICASE+10 


*  PRINT  INITIAL  DATA  * 
********************** 


C 


C 


C 


IF  (IOUT 

.EQ.  1) 

THEN 

WRITE 

(9 ,  *) 

’  Orbit  Data' 

WRITE 

(9,*) 

’  **********’ 

WRITE 

(9,*) 

WRITE 

(9.15) 

EO , AO , Y(2 , 1) 

WRITE 

(9,16) 

EF.AF, Y2FCK 

WRITE 

(9,*) 

WRITE 

(9,*) 

WRITE 

(9,*) 

’  Primary  Body  Data’ 

WRITE 

(9,*) 

’  it**************** ' 

WRITE 

(9.*)  ' 

WRITE 

(9,17) 

MU ,  R 

WRITE 

(9.*) 

WRITE 

(9,*) 

WRITE 

(9.*) 

’  Transfer  Vehicle  Data 

WRITE 

(9 ,  *) 

’  ********************* 

WRITE 

(9.*) 

WRITE 

(9,18) 

MO, ISP 

WRITE 

(9,19) 

MDOT ,  MDOTBAR 

WRITE 

(9,*) 

WRITE 

(9 ,  *) 

END  IF 

IF  (ITYP 

’.EQ.  1) 

THEN 

WRITE 

(*,*) 

’  ********************* 

WRITE 

(*,20) 

Y2F 

WRITE 

(»,*) 

’  ********************* 

WRITE 

(9,*) 

’  ********************* 

WRITE 

(9,20) 

Y2F 

WRITE 

(9,*) 

’  ********************* 

END  IF 

WRITE  («,*) 
WRITE  (*,21) 
WRITE  (*,*) 


3  2’ 
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C 


IF  (IODT.QE. 1)  THEM 
WRITE  (9,*) 

WRITE  (9,21) 

WRITE  (9 , ») 

EMDIF 


11  FORMAT 

12  FORMAT 

13  FORMAT 
+ 

14  FORMAT 

15  FORMAT 
+ 

16  FORMAT 
+ 

17  FORMAT 
♦ 

18  FORMAT 
+ 

19  FORMAT 
+ 

20  FORMAT 

21  FORMAT 


(3X, ’Vehicle  Initial  Mass:  \F9.2,’  kg’) 

(3X, 'Propellant  Mass  Flow  Rate:  ’.E13.7,'  kg/sec’) 

(3X, ’ Nondimensional  Propellant  Mass  Flow  Rate:  ’,E13.7, 

’  /sec’ ) 

(IX, 'system  Isp  currently  set  at  ’,F6.1,’  sec  ?’) 

(3X, ’Initial  Orbit  Eccentricity:  ’ ,F6. 3, 6X, ’ Initial  ’, 
’Orbit  Semimajor  Axis:  ’.F10.2,’  km  ( ’  ,E13 . 7 , ’ ) ’ ) 

(3X, 'Final  Orbit  Eccentricity:  ’ .F6.3.6X, ’Final  ’, 

’Orbit  Semimajor  Axis:  ’.F10.2,’  km  (’ .E13.7, ’) ’) 

(3X, ’Gravitational  Parameter:  ’,E13.6,’  km  /sec  ’ ,6X, 
'Radius:  ’.F9.3,’  km’) 

(3X, ’Total  Initial  Mass:  \F9.2,’  kg ’ ,6X, ’ Specif ic  ’, 
’Impulse:  ’.F9.2,’  sec’) 

(3X, ’Propellant  Mass  Flow  Rate:  ’ ,E13. 7, ’kg/sec  (’, 

E13 . 7 , ’  /sec)') 

(3X, 'AF/AO  =  ’ , E13 . 7) 

(10X, ’V-bar’ ,21X, ’e’2 OX. ’a-bar’ 18X, ’ra-bar ’ , 17X, ’rp-bar’) 


*  BEGIN  ITERATIONS  OF  STATE  VECTOR  * 
#«**#*#**#**#**#*#***»#*#*****#***«* 


VBAR=0 . DO 
25  NXT=0 
ICNT=0 

CALL  HAMING(NXT) 

IF  (NXT.EQ.O)  THEN 

WRITE(» , *)  ’HAMING  DID  NOT  INITIALIZE' 

STOP 
END  IF 

WRITE  (* , 130)  VBAB ,Y(1,1) ,y(2,l) , RABAR ( 1 ) , RPBAR ( 1 ) 
IF  (IOUT.GE. 1)  WRITE  (9,130)  VBAR,Y(1 ,1) ,Y(2, 1) , 

♦  RABAR ( 1 ) , RPBAR ( 1 ) 

50  DO  100  1=1 , IMAX 

CALL  HAMING (NXT) 

IF  (ICASE.EQ.il)  D I FF=RABARF- RABAR ( NXT ) 

IF  (ICASE.EQ.21)  DIFF=RPBAR (NXT) -RPBARF 

IF  (DABS (D IFF) .LT. 1 .D-13)  GO  TO  80 
IF  (DIFF.GT.O.DO)  GOTO  90 

WRITE  («,«) 
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C  WRITE  (« .  130)  VBAR.Yd ,NXT) , Y(2,NXT) .RABAR(NXT) , RPBAR(NXT) 

IF  (IOUT.EQ. 1)  THEN 
WRITE  (9,*) 

WRITE  (9,130)  VBAR.Yd ,NXT) ,Y(2,NXT) ,RABAR(NXT) .RPBAR(NXT) 
ENDIF 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 


C 


C 


CALL  HAMING(NXT) 

WRITE  (*,#) 

WRITE  (*  ,  130)  VBAR,Y(1,NXT) ,Y(2,NXT) .RABAR(NXT) .RPBAR(NXT) 
WRITE  (*,*) 

IF  (IOOT.EQ. 1)  THEN 
WRITE  (9 ,  *) 

WRITE  (9,130)  VBAR,Y(1 ,NXT) ,Y(2,NXT) .RABAR(NXT) .RPBAR(NXT) 
WRITE  (9,») 

ENDIF 

NXT1=NXT 

VBAR1=VBAR 

DO  55  J* 1 , 4 
Jl=5-J 

IF  (ICASE.EQ.il)  X(J1)=RABAR(NXT1) 

IF  (ICASE.EQ.21)  X(J1) =RPBAR(NXT1) 

VAR1 (Jl) =VBAR1 
VAR2 ( Jl) =Y(1 ,NXT1) 

VAR3 (J1)=Y(2, NXT 1 ) 

IF  (ICASE.EQ.il)  VAR4 ( Jl) =RPBAR(NXT1 ) 

IF  (ICASE.EQ.21)  VAR4 (Jl) =RABAR(NXT1) 

VBAR1=VBAR1-H 
NXT1=NXT1-1 
IF  (NXT1.EQ.0)  NXT 1=4 
CONTINUE 

IF  (ICASE.EQ.il)  XBAR=RABARF 
IF  (ICASE.EQ.21)  XBAR=RPBARF 
CALL  INTERP3  (XBAR.X.VARl ,P) 

VBAR=P 

CALL  INTERP3  (XBAR.X.VAR2 ,P) 

Yd  ,NXT)  =P 

CALL  INTERP3  (XBAR , X , VAR3 ,P) 

Y(2,NXT)=P 

CALL  INTERP3  (XBAR.X, VAR4 ,P) 

IF  (ICASE.EQ.il)  RPBAR (NXT) =P 
IF  (ICASE.EQ.21)  RABAR(NXT) =P 

IF  (ICASE.EQ.il)  RABAR(NXT) =RABARF 
IF  (ICASE.EQ.21)  RPBAR(NXT) =RPBARF 

WRITE  (»  ,  130)  VBAR.Yd , NXT) ,Y(2, NXT) .RABAR(NXT) .RPBAR(NXT) 
WRITE  (»,») 
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IF  (I00T.GE.1)  THEN 

WRITE  (9,130)  VBAR.Yd ,NXT) ,Y(2,NXT) .RABAR(NXT) .RPBAR(NXT) 
WRITE  (9,*) 

END  IF 

GO  TO  150 


90  IF  (IPNT.EQ.O)  GO  TO  100 

IF  ((I/IPNT»IPNT) .EQ.I)  THEN 

WRITE  (» , 130)  VBAR,Y(1,NXT) ,Y(2,NXT) .RABAR(NXT) .RPBAR(NXT) 
IF  (IOUT.EQ. 1)  WRITE  (9,130)  VBAR,Y(1 ,NXT) ,Y(2,NXT) , 

+  RABAR(NXT) , RPBAR(NXT) 

ENDIF 


100  CONTINUE 


WRITE  (*,*)  ’Maximum  number  of  iterations  reached.’ 
WRITE  (»,*)  'Program  Terminated’ 

GOTO  170 


130  FORMAT  (3X.E20 . 13 , 3X.E20 . 13 ,3X ,E20 . 13 .3X.E20 . 13 , 3X.E20 . 13) 
135  FORMAT  (3X, ’de/dVbar  =  ’.E20.13) 

150  CONTINUE 


*  TRANSFER  TOTAL  DELTA-V  AND  TIME  CALCULATION  * 
*#*#**##**#»#*##»»#*##*##*########****##**##*## 

DELV=VBAR*A0*MAN0M0 
DEL VM=DELV* 1000. 

TOF= 1 . /MDOTBAR* ( 1 . -DEXP (-DELV/ (9.8  ID-3* ISP) ) ) /3600 . /24 . 
T0FY=T0F/365 . 25 


*  Comparison  Data  * 


•  Spiral  Transfer 

VCS1=DSQRT (MU/AO) 

VCS2=DSQRT (MU/A0/Y2F) 

DELVSP=VCS 1 - VCS2 
DELVSM=DELVSP# 1000 . 

VSPBAR=DELVSP/AO/MANOMO 

T0FSP= 1 . /MDOTBAR* ( 1 . -DEXP ( -DELVSP/ (9 . 8  ID-3* ISP) ) ) /3600 . /24 . 
T0FSPY=T0FSP/365 . 25 
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*  Hohmann  Transfer 


V1=DSQRT (2»MU* ( 1 . /AO- 1 . / (A0+A0«Y2F) ) ) 
V2=DSQRT (2*MU* ( 1 . /A0/Y2F-1 . / (A0+A0*Y2F) ) ) 
DELVHM= (VI -VCS 1 ) ♦ (VCS2-V2) 

DELVHMM=DEL VHM* 1000 . 

VHMB AR= DEL VHM/ AO / MANOMO 


WRITE  («,*) 

WRITE  (* .  160)  DELV , DELVM , VBAR 
WRITE  (*,165)  TOF.TOFY 
WRITE  (»,*) 

WRITE  (*,») 

WRITE  (*,*)  ’  Comparison  Transfer  Data’ 

WRITE  (*,*)  '  »******#****»****»**##»#’ 

WRITE  (*,») 

WRITE  (*,*)  ’  Spiral  Transfer* 

WRITE  (*,*) 

WRITE  (*,160)  DELVSP , DELVSM , VSPBAR 
WRITE  (*,165)  TOFSP .TOFSPY 
WRITE  (»,*) 

WRITE  (*,»)  ’  Hohmann  Transfer' 

WRITE  (*,») 

WRITE  (*,160)  DEL VHM , DELVHMM , VHMB AR 
WRITE  (»,») 

WRITE  (»,») 

WRITE  (*,*) 

IF  (IOUT.GE. 1)  THEN 
WRITE  (9,*) 

WRITE  (9,160)  DELV, DELVM, VBAR 
WRITE  (9,165)  TOF.TOFY 
WRITE  (9,*) 

WRITE  (9,*) 

WRITE  (9,*)  *  Comparison  Transfer  Data’ 

WRITE  (9 , »)  ’  **************»**#******’ 

WRITE  (9 , ») 

WRITE  (9 , »)  ’  Spiral  Transfer’ 

WRITE  (9 , ») 

WRITE  19,160)  DELVSP .DELVSM, VSPBAR 
WRITE  (9,165)  TOFSP, TOFSPY 
WRITE  (9.*) 

WRITE  (9,*)  ’  Hohmann  Transfer’ 

WRITE  (9,*) 

WRITE  (9,160)  DELVHM, DELVHMM, VHMB AR 
WRITE  (9,*) 

WRITE  (9,*) 

WRITE  (9,*) 

END  IF 
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IF  (ITYP.EQ. 1)  THEN 
WRITE  (10,*)  Y2F.VBAR 
WRITE  (11.*)  Y2F.VSPBAR 
WRITE  (12,*)  Y2F , VHMBAR 
END  IF 
C 

160  FORMAT  (3X,’Total  Transfer  Delta-V:  ’.F8.4,’  km/s  (’.F8.2, 
+  ’  m/s) ’ ,6X, ’Nondimensional  Delta-V:  ’,F8.4,’  /sec’) 

165  FORMAT  (3X, ’Total  Transfer  Time:  \F9.3,’  days  (’, 

*  F7 . 3 , '  yr) ’ ) 


DIFFYF=Y2FCK-Y2F 

IF  (DABS(DIFFYF) .GT. l.E-13)  GOTO  10 


170  CONTINUE 


IF  (IOUT.GE. 1)  CLOSE  (9) 
IF  (ITYP.EQ. 1)  THEN 
CLOSE  (10) 

CLOSE  (11) 

CLOSE  (12) 

END  IF 


STOP 

END 


********************************* 
*  SUBROUTINE  RHS  * 

»»  i  «  *  »  i  **«»**»**»***»***«**»»**»  » 

SUBROUTINE  RHS(NXT) 

IMPLICIT  DOUBLE  PRECISION  (A-H.0-Z) 

COMMON  /HAM/  VBAR, Y(42 , 4) ,F (42 , 4) .ERR (42) ,N .H.MODE 
COMMON  /RH/  RABAR(4) ,RPBAR(4) , ICASE , ICNT , DIFF 


CALL  I NTERPOA ( Y ( 1 , NXT ) , DELAA , DELEA , ICNT) 
CALL  INTERPOP (Y( 1 , NXT) , DELAP , DELEP , ICNT) 
ICNT=ICNT+ 1 

*  Define  Derivatives 

PI=DBLE (ACOS (-1.0)) 

•  C 

IF  (ICASE.LT. 20)  THEN 
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DELE=DELEP 
DEL A= DELAP 

F ( 1 ,NXT) =Y (2 ,NXT) ** . 5« (DELEP-DELEA) / (4 . »PI) 

F (2 ,  NXT) =Y(2 , NXT) **1.5* (DELAP-DELAA) / (4 . *PI ) 

ENDIF 

C 

IF  (ICASE.GT. 20)  THEN 
DELE=DELEA 
DEL A- DEL AA 

F { 1 ,NXT) aY (2 ,NXT) »* .  5* (DELEA-DELEP) / (4 . *PI) 

F (2 ,NXT) =Y (2 ,NXT) »* 1 . 5* (DELAA-DALAP) / (4 . *PI ) 

END  IF 
C 

IF  (F ( 1 .NXT) .LT.0.D0)  THEN 
F ( 1 ,NXT)=Y(2,NXT)*».5*DELE/(2.*PI) 

F (2 ,NXT) =Y(2 ,NXT) #*  1 . 5#DELA/ (2 . *PI ) 

WRITE  (*,50)  Y ( 1 , NXT ) , Y ( 2 , NXT ) 

END  IF 

*  Define  Changes  in  ra-bar  and  rp-bar 

RABAR (NXT) =Y (2 ,  NXT) *  ( 1 . +Y( 1 ,NXT) ) 

RPBAR(NXT) =Y(2 , NXT) *  ( 1 . -Y( 1 , NXT) ) 

C 

50  FORMAT  (3X,’ Single  Orbit  Performed  at  e  =  ’,E13.7,’  and  abar 
♦  E13.7) 

RETURN 

END 

C 

C 

SUBROUTINE  HAMING(NXT) 

C  VERSION  OF  11/20/1987 
C  PURPOSE 

C  HAMING  IS  AN  ORDINARY  DIFFERENTIAL  EQUATIONS  INTEGRATOR 
C  IT  IS  A  FOURTH  ORDER  PREDICTOR-CORRECTOR  ALGORITHM  WHICH 
C  MEANS  THAT  IT  CARRIES  THE  LAST  FOUR  VALUES  OF  THE  STATE 

C  VECTOR,  AND  EXTRAPOLATES  THESE  VALUES  TO  OBTAIN  A  PREDICTED 

C  NEXT  VALUE  (THE  PREDICTION  STEP)  AND  EVALUATES  THE  EQUATIONS 
C  OF  MOTION  AT  THE  PREDICTED  POINT,  AND  THEN  CORRECTS  THE 

C  EXTRAPOLATED  POINT  USING  A  HIGHER  ORDER  POLYNOMIAL  (THE 

C  CORRECTION  STEP) . 

C  INPUT 

C  NXT  --  IN  THE  CALL  SPECIFIES  WHICH  OF  THE  FOUR  VALUES  OF 
C  THE  STATE  VECTOR  IS  THE  CURRENT  ONE.  NXT  IS  UPDATED 

C  BY  HAMING  AUTOMATICALLY,  BUT  MUST  BE  SET  TO  ZERO  ON 

C  THE  FIRST  CALL. 

C  CALL  ROUTINES 
C  RHS(NXT) 

C  REFERENCES 
C  WILLIAM  WEISEL 
C  PROGRAMMER 
C  RODNEY  D.  BAIN 
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C  PROGRAM  MODIFICATIONS 
C  NONE 
C  COMMENTS 

C  TOL  —  IS  HAMING'S  START  DP  TOLERANCE  ...  SET  TO  REASONABLE  VALUE 

C  AS  NECESSARY 

C  THE  COMMON  BLOCK  CONTAINS: 

C  VBAR  —  IS  THE  INDEPENDENT  VARIABLE  (OFTEN  TIME) 

C  Y(42 ,4)  —  IS  THE  STATE  VECTOR,  4  COPIES  OF  IT,  WITH  NXT  POINTING 

C  POINTING  TO  THE  CURRENT  ONE,  THE  LIMIT  OF  42  EQUATIONS 

C  OF  MOTION  CAN  BE  CHANGED. 

C  F (42, 4)  —  ARE  THE  EQUATIONS  OF  MOTION  EVALUATED  AT  THE  SAME  TIMES 
C  AS  THE  STATE  VECTOR  Y  ...  IT  IS  THE  JOB  OF  SUBROUTINE 

C  RHS  TO  CALCULATE  THESE. 

C  ERR (42)  --  IS  AN  ESTIMATE  OF  THE  ONE-STEP  INTEGRATION  ERROR 

C  N  --  IS  THE  NUMBER  OF  ODES  . . .  LIMIT  IS  42  UNLESS  YOU  CHANGE 

C  THE  COMMON  BLOCK 

C  H  --IS  THE  TIMESTEP  ...  ONE  CALL  TO  HAMING  INCREMENTS  X  BY  H 

C  MODE  --  IS  ZERO  FOR  EOM  ONLY,  1  FOR  EOM  AND  EOV 

C  THE  USER  MUST  SUPPLY  A  MAIN  PROGRAM,  AND  THE  SUBROUTINE  RHS (NXT)  WHICH 
C  EVALUATES  THE  EQUATIONS  OF  MOTION. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  /HAM/  VBAR,Y(42 ,4) ,F(42 ,4) ,ERR(42) ,N,H,MODE 
DATA  ZERO , ONE , TWO , THREE , FOUR/O . DO , I . DO , 2 . DO , 3 . DO , 4 . DO/ 

TOL= 1 . D- 12 

CHECK  IF  THIS  IS  THE  FIRST  CALL  . . .  HAMING  (LIKE  ALL  PREDICTOR- 
CORRECTORS)  NEEDS  ’PREVIOUS'  VALUES 

IF (NXT)  190,10,200 
C 

C  IT  IS  A  PICARD  INTERATION  (SLOW  AND  EXPENSIVE)  TO  STEP  BACKWARDS 

C  IN  TIME  THREE  STEPS  TO  GET  THE  4  PREVIOUS  POINTS.  A  SUCCESSFUL 

C  STARTUP  RETURNS  NXT=1,  AND  TIME  HAS  NOT  BEEN  INCREMENTED.  IF 

C  STARTUP  FAILS,  NXT  WILL  BE  RETURNED  AS  ZERO. 

C 

10  XO=VBAR 
HH=H/TWO 
CALL  RHS ( 1 ) 

DO  40  L=2 ,4 
VBAR=VBAR+HH 
DO  20  1*1, N 

Y(I,L)=Y(I,L-1)+HH«F(I,L-1) 

20  CONTINUE 

CALL  RHS (L) 

VBAR=VBAR+HH 
DO  30  1=1, N 

Y(I ,L) =Y( I ,L- 1) +H*F ( I ,L) 

30  CONTINUE 

CALL  RHS (L) 

40  CONTINUE 
JSW=-10 


o  o  o  o 


50  ISW=1 

DO  120  1  =  1, N 

HH=Y(I , 1) +H* (9 .D0*F (I , 1 ) + 19 . D0»F (I , 2) -5 . D0»F (I , 3) 

1  +F(I,4))/24.D0 

IF(DABS(HH-Y(I ,2) ) .LT.TOL)  GOTO  70 
ISW=0 

70  Y(I , 2) =HH 

HH=Y(I , 1) +H» (F (I , 1) +FOUR«F (I , 2) +F (I ,3) ) /THREE 

IF (DABS (HH-Y( I, 3)). LT.TOL)  GOTO  90 

ISW=0 

90  Y(I ,3) =HH 

HH=Y(I , 1) +H# (THREE»F (I , 1) +9 . DO*F (I , 2) +9 . DO«F ( I ,3) 

1  +THREE*F (I , 4) ) /8 . DO 
IF(DABS(HH-Y(I .4)) .LT.TOL)  GOTO  110 
ISW=0 

110  Y(I ,4) =HH 
120  CONTINUE 
VBAR=XO 
DO  130  L=2 , 4 
VBAR=VBAR+H 
CALL  BHS(L) 

130  CONTINUE 

IF(ISW)  140,140,150 
140  JSW=JSW+1 

IF(JSW)  50,280,280 
150  VBAR=XO 
ISW=1 
JSW=1 

DO  160  1  =  1, N 
ERR(I) =ZERO 
160  CONTINUE 
NXT=  1 
GOTO  280 
C 

C  A  CALL  TO  HAMING  WITH  NXT=-NXT,  AFTER  A  SUCCESSFUL  STARTUP, 

C  WILL  TURN  OFF  THE  SECOND  EVALUATION  OF  THE  EQUATIONS  OF  MOTION 
C  FOLLOWING  THE  CORRECTOR  STEP.  IN  SYSTEMS  WHERE  THE  EQUATIONS  OF 
C  MOTION  ARE  VERY  EXPENSIVE,  THIS  CAN  HALVE  YOUR  RUN  TIME. 

C 

190  JSW=2 

NXT=IABS (NXT) 

THIS  IS  THE  PREDICTOR-CORRECTOR  ALGORITHM  . . .  FIRST  THE  INDICES 
ARE  PERMUTED 

200  VBAR=VBAR+H 

NP1=M0D(NXT,4)+1 
GOTO  (210,230) ,ISW 
210  GOTO  (270,270,270,220) , NXT 
220  ISW=2 

230  NM2=M0D (NP1 , 4) + 1 
NM1=M0D(NM2 ,4) +1 
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NP0=M0D(NM1 , 4) +1 


. . .  THEN  THE  PREDICTOR  PART  IS  RUN  TO  FIND  AN  EXTRAPOLATED  VALUE 
OF  THE  STATE  VECTOR  AT  THE  NEW  TIME  ... 

DO  240  1=1, N 

F (I ,NM2) =Y(I ,NP1) +F0UR#H» (TW0*F(I ,NP0) -F(I ,NM1 ) 

1  +TW0»F (I ,NM2) ) /THREE 

Y(I,NP1)=F(I,NM2)-0.925619835D0*ERR(I) 

240  CONTINUE 

THE  EQUATIONS  OF  MOTION  ARE  EVALUATED  AT  THE  EXTRAPOLATED  VALUE 
OF  THE  STATE  VECTOR  . . . 

CALL  RHS(NPl) 

AND  THE  CORRECTOR  ALGORITHM  IS  USED  TO  ADD  THIS  NEW  INFORMATION 
AND  OBTAIN  A  BETTER  VALUE  OF  THE  NEW  STATE  VECTOR  . . . 


DO  250  1=1, N 

Yd  ,NP1 )  =  (9 . D0*Y(I , NPO)  —  Y ( I ,  NM2)  +THREE«H*  (F (I  ,NP1 ) 

1  +TWO#F (I ,NP0) — F ( I , NM1) ) ) /8 .DO 

ERR ( I ) =F ( I , NM2 ) - Y ( I , NP 1 ) 

Y(I,NPl)=Y(I,NPl)+0. 0743801 653D0»ERR( I) 

250  CONTINUE 

GOTO  (260,270) ,JSW 

FINALLY,  THE  EQUATIONS  OF  MOTION  ARE  RE-EVALUATED  AT  THE  BETTER 
VALUE  OF  THE  STATE  VECTOR  ...  THIS  CAN  BE  SUPPRESSED. 


260  CALL  RHS(NPl) 

270  NXT=NP1 
280  RETURN 
END 
C 
C 

C  ##*#######*###**#*#*##**»######## 
C  «  SUBROUTINE  INTERP3  * 

C  *  * 

C  *  The  following  subroutine  interpolates  between  four  data  * 
C  *  points  using  a  third  order  polynomial.  * 

C  ***»»»»»*»»»*»»*»«****««»««««»»»* 

c 

SUBROUTINE  INTERP3  (XBAR,X,D,P) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  X(4) ,P1 (3) ,D(4) 

C 

DO  300  K=1 ,3 

DO  300  J= (K+l)  ,4 
Jl=4-J+ (K+ 1 ) 

D(J1)=(D(J1)-D(J1-1))/(X(J1)-X(J1-K)) 
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300  CONTINUE 
C 

Z=XBAR-X(1) 

PI  ( 1) *D( 1) +D(2) *Z 
DO  310  J=2 ,3 
Z=Z* (XBAR-X(J) ) 

PI (J) *P1 (J-l) +D(J+1) »Z 
310  CONTINUE 
C 

P=P1 (3) 

C 

RETURN 

END 

C 

C 

C 

C 

C  »##*#***#########*»#*#####*#«## 

C  *  SUBROUTINE  INTERPO  * 

C  »##»#*###»###*»*#*#*####*»**#*# 

c 

c  »##*#*###**###»***####*##**#*#* 

C  »  The  following  subroutine  utilizes  the  Newton  formula  * 

C  *  to  interpolate/extrapolate  from  the  given  data  sets  * 

C  *  utilizing  an  4th  degree  interpolating  polynomial.  * 

C  a##*#*#**#****#*######*#####*## 

C 
C 

SUBROUTINE  INTERPOA (EBAR , DA , DE , ICNT) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E (100) ,LAMDA(100) ,DELA(100) ,DELE(100) 

C 

C 

C  *  OPEN  DATA  FILES  * 

C  »»»*»»**»»*»»»»»*»» 

c 

C  *  Constant  ra  Data 

C 

IF  (ICNT.EQ.O)  THEN  ,,,, 

OPEN  (UNIT* 15,  STATUS  =  ’OLD’ .FILE  =  ’ALAM.DAT’) 

OPEN  (UNIT* 16 ,  STATUS  *  ’OLD’, FILE  *  ’ADELA.DAT’) 

OPEN  (UNIT* 17,  STATUS  *  ’OLD’ .FILE  =  ’ADELE.DAT’) 

I DATA* 99 
END  IF 


C 


IF  (ICNT.EQ.O)  THEN 

REWIND  (15) 

REWIND  (16) 

REWIND  (17) 
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»  INPUT  INITIAL  DATA  * 
a*#####**##*######**## 


DO  5,  Is 1 .IDATA 

READ  (15, «}  E(I) ,LAMDA(I) 
READ  (16, *)  E(I),DELA(I) 
READ  (17, «)  E(I) , DELE (I) 

5  CONTINUE 

IDEG=4 

CLOSE  (15) 

CLOSE  (16) 

CLOSE  (17) 

ENDIF 


*  EVALUATE  FUNCTION  AT  SPECIFIED  VALUES  OF  ECCENTRICITY  » 
a**##*######*#*#*#**####*****####***#####*#**######**##** 


*  Check  location  of  data  point  * 

DO  30  1=1, IDATA 

DIFF=EBAR-E(I) 

IF  (ABS(DIFF) .LE.l.E-15)  THEN 
LAM=LAMDA(I) 

DA=DELA(I) 

DE=DELE (I) 

GO  TO  40 
ENDIF 

IF  ((I.EQ.l) .AND. (DIFF.LT.O.))  THEN 
IMIN= 1 

IMAX=IMIN+IDEG 
GO  TO  35 
ENDIF 

IF  ( (I. EQ. IDATA) .AND. (DIFF.GT.O. ) )  THEN 
I  MAXIDATA 
I MI N= I MAX- I DEG 
GO  TO  35 
ENDIF 

ILEFT=I 

IF  (DIFF.LT.O.)  GO  TO  31 
30  CONTINUE 
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31  IDEGP1=IDEG+ 1 

IDEGL=IDEGPl/2 
I MI N* I LEFT+ 1 - I DEGL 

•  I  MAX55 1  MI  N+ 1  DEG 
C 

IF  (IMIN.LT. 1)  THEN 
IMIN= 1 

IMAX= IMIN+IDEG 
END  IF 

•  C 

IF  ( I MAX. GT. I DATA)  THEN 
I  MAXIDATA 
I MI N= I MAX- I DEG 
END  IF 


*  Evaluate  Polynomial  * 

35  CALL  INTERP(IDEG,IMIN,IMAX,EBAR,E,LAMDA,DELA,DELE,LAM,DA,DE) 

40  CONTINUE 
RETURN 
END 


SUBROUTINE  INTERPOP (EBAR . DA . DE , ICNT) 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E (100) .LAMDA(IOO) .DELA(IOO) .DELE(IOO) 


#  OPEN  DATA  FILES  » 
»»»*»*»»»»»*»»»»»»» 


«  Constant  rp  Data 

IF  (ICNT.EQ.O)  THEN 

OPEN  (UNIT= 18 ,  STATUS  = 
OPEN  (UNIT55 19 ,  STATUS  = 
OPEN  (UNIT55 20  ,  STATUS  = 
IDATA=98 
END  IF 


IF  (ICNT.EQ.O)  THEN 

REWIND  (18) 

REWIND  (19) 

REWIND  (20) 


’OLD’,  FILE  55  ’PLAM.DAT’ ) 
’OLD’, FILE  =  'PDELA.DAT’) 
'OLD '.FILE  55  'PDELE.DAT') 
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*  INPUT  INITIAL  DATA  # 
####»##»###*#*######*# 


C 

C 


c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


c 


c 

c 


DO  5,  1=1, IDATA 

READ  (18, *)  E(I) ,LAMDA(I) 
READ  (19,#)  E(I),DELA(I) 
READ  (20,#)  E ( I ) , DELE (I) 

5  CONTINUE 

IDEQ=4 

CLOSE  (18) 

CLOSE  (19) 

CLOSE  (20) 

ENDIF 


»  EVALUATE  FUNCTION  AT  SPECIFIED  VALUES  OF  ECCENTRICITY  » 
*»#***#*****##*»##*###**»#**#*#****»##****#*#*##******##* 


#  Check  location  of  data  point  # 


DO  30  1=1, IDATA 

DIFF=EBAR-E(I) 

IF  (ABS(DIFF) .LE. l.E-15)  THEN 
LAM=LAMDA(I) 

DA=DELA(I ) 

DE=DELE(I) 

GO  TO  40 
ENDIF 

IF  ((I.EQ. 1) .AND. (DIFF.LT.O.))  THEN 
IMIN= 1 

IMAX=IMIN+IDEG 
GO  TO  35 
ENDIF 

IF  ((I.EQ. IDATA) .AND. (DIFF.GT.O. ) )  THEN 
IMAX=IDATA 
IMIN=IMAX-IDEG 
GO  TO  35 
ENDIF 

ILEFT=I 

IF  (DIFF.LT.O.)  GO  TO  31 
30  CONTINUE 


31  IDEGP1=IDEG+ 1 
IDEGL=IDEGPl/2 
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o  o  o  o 


I MI N= I LEFT  + 1 - 1 DEGL 

IMAX=IMIN+IDEG 

C 

•  IF  (IMIN.LT. 1)  THEN 

IMIN= 1 

IMAX=IMIN+IDEG 
END  IF 
C 

IF  (IMAX.GT. IDATA)  THEN 

•  I MAX3 IDATA 

IMIN=IMAX-IDEG 
END  IF 


*  Evaluate  Polynomial  * 


C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 


C 


35  CALL  INTERP ( IDEG , IMIN , I MAX , EBAR , E , LAMDA , DELA , DELE , LAM, DA , DE) 

40  CONTINUE 
RETURN 
END 


a*###*#*###*#*#*##*#**##*#***#*# 

»  SUBROUTINE  INTERP  » 

*##»»##*####*#*»###»#####*##**## 

a###################*#*########* 

•  The  following  algorithm  evaluates  an  interpolating  * 

«  nth  degree  polynomial  « 

»#*#*##*#*#»*#»##**#*###*###»### 

SUBROUTINE  I NTERP ( I DEG , I MI N , I MAX , EB AR , E , L AMD A , DEL A , DELE , L AM , DA , DE ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.L-Z) 

DIMENSION  E ( 100)  , LAMDA ( 100)  .DELAUOO)  ,DELE(100) 

DIMENSION  X(IOO) ,D(100) ,P(100) 

IDEGP1=IDEG+1 

DO  200  1=1,3 

DO  150  J=IMIN, IMAX 
K=J-IMIN+ 1 

IF  (I.EQ.l)  D(K) =LAMDA(J) 

IF  (I.EQ.2)  D (K) =DELA(J) 

IF  (I.EQ.3)  D (K) =DELE (J) 

X(K)=E(J) 

150  CONTINUE 

DO  155  K= 1 , IDEG 
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DO  155  J= (K+l) , IDEQP1 
J1=IDEGP1-J+ (K+ 1) 

D(«J1)  =  (D(J1)  — D(J1  —  1) )  /  (XCJ1)  -X(Jl-K) ) 
CONTINUE 

Z=EBAR-X(1) 

P  ( 1)  =D  ( 1)  +D  (2)  *Z 

DO  160  J=2 , IDEG 
Z=Z* (EBAR-X (J) ) 

P(J)=P(J-1)+D(J+1)»Z 

CONTINUE 

IF  (I.EQ.l)  LAM=P (IDEG) 

IF  (I.EQ.2)  DA=P(IDEG) 

IF  (I.EQ.3)  DE=P(IDEG) 

CONTINUE 

RETURN 

END 
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Appendix  D-l: 


Program  TRANSMUL  Data  -  LEO  to  GEO  Transfer 


Orbit  Data 


Initial  Orbit  Eccentricity:  0.000 
Final  Orbit  Eccentricity:  0.000 


Initial  Orbit  Seaiaajor  Axis:  6678.14  ka  (0. 1000000E+01) 

Final  Orbit  Seaiaajor  Axis:  42241.15  ka  (0.6325281E*01) 


Priaary  Body  Data 
Gravitational  Paraaeter: 


3  2 

398,601  ka  /see 


Badiue:  6378.145  ka 


Trane fer  Vehicle  Data 

Total  Initial  Maea:  5000.00  kg  Specific  Impulse:  863.00  eec 

Propellant  Haas  Floe  Bate:  0.3977100E-03kg/aec  (0.7954200E-07  /sec) 


Bey 

Tiae  (aec) 

IU  (Bad! 

V-bar 

e 

a-bar 

ra-bar 

rp-bar 

1 

O.OOOOE+OO 

0.000 

O.OOOOE+OO 

0.0000E+00 

0. 1000E+01 

0.1000E+01 

0. 1000E+01 

50 

0.2723E+06 

2.196 

0.2400E-01 

0.2568E-01 

0.1026E+01 

0.1052E+01 

0.1000E+01 

100 

0 . 5607E+06 

3.591 

0.5000E-01 

0.5348E-01 

0.1056E+01 

0.1113E+01 

0, 1000E+01 

150 

0.8637E+06 

3.084 

0.7800E-01 

0.8333E-01 

0.1090E+01 

0.1181E+01 

0.1000E+01 

201 

0.1190E+07 

0.554 

0.1090E+00 

0.1162S+00 

0.1131E+01 

0.1263E+01 

0.1000E+01 

251 

0.1527E+07 

1.316 

0.142QE+Q0 

0.1511E+00 

0. 1178E+01 

0.1356E+01 

0.1000E+0! 

301 

0. 1894E+07 

0.507 

0. 1790E+00 

0.1899E+00 

0.1234E+01 

0.1460E+01 

0.1000E+01 

349 

0.2277E+07 

8.258 

0 . 2 190E+00 

0.2315E+00 

0. 1301E+01 

0.1602E+01 

0.1000E+01 

400 

0.2655E+07 

4.517 

0.2800E+00 

0.2736E+00 

0.1376E+01 

0. 1753E+01 

0.1000E+01 

450 

0 . 3054E+07 

4.556 

0.3050E+00 

0.3190E+00 

0. 1468E+01 

0.1937E+01 

0.1000E+01 

501 

0.3420E+07 

0.235 

0.348OE+OO 

0.3817E+00 

0. 1568E+01 

0.2133E+01 

0. 1000E+01 

550 

0 . 3828E+07 

0.535 

0 . 3980E+00 

0.4102E+00 

0. 1695E+01 

0.2391E+01 

0.1000E+01 

600 

0.4218E+07 

2.889 

0.4480E+00 

0.4574E+00 

0. 1843E+01 

0.2686E+01 

0. 1000E+01 

650 

0.4813E+07 

0.827 

0.5010E+00 

0.5058E+00 

0.2023E+01 

0.3047E+01 

0.1000E+01 

700 

0.4995E+07 

1.249 

0.5550E+Q0 

0.5533E+00 

0.2238E+01 

0.3477E+01 

0.1000E+01 

750 

0 . 5432E+07 

1.126 

0.6200E+00 

0.8078E+00 

0.2550E+01 

0.4100E+01 

0.1000E+01 

775 

0.5669E+07 

3.036 

0. 6570E+00 

0.8375E+00 

0.2758E+01 

0.4517E+01 

0.1000E+01 

800 

0.5892E+07 

8.011 

0 . 6930E+00 

0.6853E+00 

0.2988E+01 

0.4977E+01 

0.1000E+01 

825 

0 . 6102E*07 

3.337 

0.7280E+Q0 

0.6914E+00 

0.3241E+01 

0.5482E+01 

0.1000E+01 

849 

0 . 6368E+07 

5.017 

0.7740E+00 

0.7243E+00 

0.3627E+01 

0.6254E+01 

0. 1000E+01 

850 

0.6389E+07 

4.715 

0.8085E+00 

0 . 7269E+00 

0.3662E+01 

0.6325E+01 

0.1000E+01 

855 

0.8552E+07 

1.431 

0.8378E+00 

0.8710E+00 

0.3785E+01 

0.6325E+01 

0. 1245E+01 

860 

0 . 6765E+07 

1.006 

0.8772E+00 

0.5884E+00 

0.3982E+01 

0.6325E+01 

0. 1638E+01 

865 

0 . 6992E+07 

0.781 

0 . 9208E+00 

0.4893E+00 

0.4247E+01 

0.6325E+01 

0.2168E+01 

870 

0.7250E+07 

3.835 

0.9728E+00 

0.3624E+00 

0.4642E+01 

0.6325E+01 

0.2959E+01 

875 

0.7528E+07 

0.281 

0. 1031E+01 

0.2119E+00 

0.5219E+01 

0.6325E+01 

0.4112E+01 

879 

0 . 7878E+07 

4.915 

0.11 10E+01 

-0.3053E-10 

0.6325E+01 

0.6325E+01 

0.6325E+01 
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Total  Transfer  Delta-?:  8.5780  ka/s  (  8577.98  a/s) 

Total  Transfer  Tiae:  92.682  days  (  0.25  yr) 

Coapariaon  Transfer  Data 
Spiral  Transfer 

Total  Transfer  Delta-?:  4.6539  ka/s  (  4653.90  a/s) 

Total  Transfer  Tiae:  61.534  days  (  0.17  yr) 

Hohaann  Transfer 

ToUl  Transfer  Delta-?:  3.8937  ka/s  (  3893.75  a/s) 


Appendix  D-2:  Program  TBANSALT  Data  z.  GEO  Transfer 


Orbit  Data 


Initial  Orbit  Eccentricity: 
Final  Orbit  Eccentricity: 

0.000 

0.000 

Initial  Orbit  Seaiaajor  Axis: 
Final  Orbit  Seaiaajor  Axis: 

6678.14  ka 
42241.15  ka 

(0.10000008+01) 

(0.6325281E+01) 

Priaary  Body  Data 

0 

Gravitational  Paraaeter: 

* 

398,601  ka  /sec  Radius:  6378.145  ka 

Transfer  Vehicle  Data 


Total  Initial  Mass:  5000.00  kg  Specific  Iapulse:  863.00  sec 

Propellant  Mass  Flow  Hate:  0.39771001-03  kg/sec  (0.7954203E-07  /sec) 


V-bar 

0.0000000000E+00 
0.5000000000E-01 
0.10000000008+00 
0. 1500000000E+00 
0. 2000000000E+00 
0.2500000000E+00 
0 . 3000000000K*00 
0 . 3500000000E+00 
0.40000000008+00 
0. 4500000000E+00 
0.50000000008+00 
0 . 5500000000E+00 
0.6000000000E+00 
0. 85000000008*00 
0 . 7000000000E+00 
0.75000000008*00 
0.8000000000E*00 
0 . 8500000000E+00 
0.90000000008*00 
0.9500000000E*00 
0.1000000000E+01 
0.10500000008*01 
0.11000000008*01 
0.11250949758*01 


e 

o.ooooooooooe*oo 

0.13025314798-08 
0.2603426238E-08 
0. 39026401 18E-08 
0.5200126439E-08 
0. 649583580 IE-08 
0. 77897 158098*08 
0.90817108098-08 
0.10371761568-07 
0.11659804888-07 
0.1294577325E-07 
0.14229594298-07 
0.15511190298-07 
0.16790477538-07 
0.1806738559E-07 
0.1934175848E-07 
0.20813543588-07 
0.2188261064E-07 
0.2314883016E-07 
0.2441208177E-07 
0.25672150178-07 
0.2892892257E-07 
0.2818218550E-07 
0.2880980737E-07 


a-bar 

0. 1000000000E*01 
0.10557704068+01 
0.11163398228+01 
0.11822750688+01 
0.12542291908+01 
0.13329576908+01 
0.14193384428+01 
M514396269E+01 
0. 16193334798+01 
0.17355680738+01 
0.18647818978*01 
0.20089818298+01 
0.21705781758+01 
0.23524860468+01 
0.25582577568+01 
0.27922575728+01 
0.30598950488+01 
0.33679404938+01 
0.3724957392E+01 
0.4141904050E+01 
0.46329846278+01 
0.5216875095E+01 
0.59185254898+01 
0.6325281014E+01 


ra-bar 

0.10000000008+01 
0.10557704088+01 
0.1I18339825E+0I 
0.11822750738+01 
0.12542291968+01 
0.13329576998+01 
0.1419338453E+01 
0. 1514396283E+01 
0.16193334968+01 
0.17355680938+01 
0.18647819218+01 
0. 200898 1858E+01 
0.21705782098+01 
0.23524860858+01 
0.25582578028+01 
0.27922576268+01 
0.30598951118+01 
0.33679405668+01 
0.3724957478E+01 
0.414 1904 15 1E+0 1 
0.46329847466+01 
0.5218875236E+01 
0.5918525656E+01 
0.63252811978+01 


rp-bar 

0.1000000000E+01 
0. 1055770405E+01 
0.II16339820E+01 
0. 1182275064E+01 
0. 1 254229 183E+0 1 
0. 133295768 1E+01 
0.14193384318+01 
0.1514396255E+01 
0.16193334628+01 
0.1735568052E+01 
0.18647818738+01 
0.2008981801E+01 
0.21705781418+01 
0.23524880068+01 
0.25582577098+01 
0.27922575188+01 
0.30598949856+01 
0.33879404198+01 
0.37249573068+01 
0.41419039496+01 
0.46329845098+01 
0.5216874955E+01 
0 . 59185253228+01 
0.6325280832E+01 


Total  Transfer  Delta-7:  8.8922  ka/i  (  8892.22  ■ /•)  londiaensional  Delta-7:  1.1251  /sec 

Total  Transfer  Tiae:  93.390  days  (  0.256  yr) 

Cogarison  Transfer  Data 
Spiral  Transfer 

Total  Transfer  Delta-7:  4.6539  ka/s  (  4653.90  a/s)  londiaensional  Delta-7:  0.6024  /sec 

Total  Transfer  Tine:  61.534  days  (  0.168  yr) 

Bohaann  Transfer 

Total  Transfer  Delta-7:  3.8937  ka/s  (  3893.75  n/s)  londiaensional  Delta-7:  0.5040  /sec 
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Appendix  D-3:  Program  TRANSMUL  Data  -  Transfer  to  Molniya  Orbit 


Orbit  Data 

Initial  Orbit  Eccentricity: 
Final  Orbit  Eccentricity: 


Initial  Orbit  Seaiaajor  Axis: 
Final  Orbit  Seaiaajor  Axis: 


7184.76  ka  (0. 1000000E+01) 
26610.23  ka  (0 . 37037058*0 1 ) 


Priaary  Body  Data 

3  2 

Gravitational  Paraaeter:  398,601  ka  /sec  Radius:  6378.145  ka 


Transfer  Vehicle  Data 

Total  Initial  Mass:  5000.00  kg  Specific  Iapulse:  863.00  sec 

Propellant  Mass  Flow  Bate:  0 . 3977 100E-03kg/sec  (0.7954200E-07  /sec) 


Bev 

Tiae  (sec) 

15  LMl 

V-bar 

8 

a-bar 

ra-bar 

rp-bar 

1 

0.00008+00 

0.000 

0.0000E+0Q 

0 . OOOOE+OO 

0 . 1000E+01 

0.1000E+01 

0.1000E+01 

50 

0.3059E+06 

3.138 

0.2800E-01 

0.2996E-01 

0.10308+01 

0.1081E+01 

0. 1000E+01 

101 

0.6359E+06 

2.159 

0.5900E-01 

0. 6308E-01 

0.10678+01 

0.1134E+01 

0.1000E+01 

150 

0.9672E+06 

0.332 

0.9I00E-01 

0.9716E-0I 

0.1107E+01 

0.I215E+01 

0. 1000E+01 

200 

0 . 1338E+07 

1.554 

0. 1280E+00 

0. 1363E+00 

0. 1157E+01 

Q.1315E+01 

0.1000E+01 

250 

0. 1746E+07 

0.156 

0.1700E+00 

0 . 1805E+00 

0.12208+01 

0.1440E+01 

0.1000E+01 

300 

0.21488+07 

1.396 

0.21308+00 

0  2253E+00 

0. 1290E+01 

0.1581E+01 

0.1000E*01 

325 

0.2330E+07 

2.174 

0.2330E+00 

0 . 2459E+00 

0 . 1326E+01 

0.1652E+01 

0 . 1000E+0 1 

350 

0.25448+07 

3.604 

0. 25708+00 

0.2705E+00 

0.137QE+01 

0.1741E+01 

0.1000E+01 

376 

0.2727E+07 

0.468 

0.2780E+00 

0.2919E+00 

0.1412E+01 

0. 1824E+01 

0. 1000E+01 

400 

0.29338+07 

5.557 

0.3020E+00 

0.31808+00 

0. 14828+01 

0.19248+01 

0.10008*01 

425 

0.31348+07 

3.548 

0.3260E+00 

0.3400E+00 

0.15158+01 

0.2030E+01 

0. 1000E+01 

450 

0 . 3340E+07 

4.354 

0.35108+00 

0.3646E+00 

0.15748+01 

0.21488*01 

0. 1000E+01 

475 

0.35338+07 

0.378 

0 . 3750E+00 

0.3880E+00 

0 . 1634E+01 

0.2268E+01 

0.10008*01 

500 

0.37688+07 

5.816 

0.40508+00 

0.4169E+00 

0.1715E+01 

0.2430E+01 

0 . 1000E+01 

525 

0.4005E+07 

4.157 

0 . 4360E+00 

0.4462E+00 

0.1805E+01 

0.26118*01 

0.10008*01 

550 

0 . 4257E+07 

3.334 

0.4700E+00 

0.4777E+00 

0. 1914E+01 

0.2829E+01 

0.1000E+01 

575 

0,44748+07 

3.958 

0. 50008*00 

0.50498+00 

0.20198*01 

0.3039E+01 

0.1000E+01 

600 

0.4692E+07 

4.796 

0 . 5310E+00 

0 . 5324E+00 

0.2138E+01 

0.3277E+01 

0.1000E+01 

625 

0.4971E+07 

2.469 

0.5720E+00 

0.5678E+0 0 

0.2314E+01 

0.3628E+01 

0. IOO0E+O1 

653 

0.52478+07 

4.174 

0.81408+00 

0.60298+00 

0.2518E+01 

0.4037E+01 

0. 1000E+01 

675 

0.54898+07 

2.338 

0 . 6490E+00 

0.631 1E+00 

0.2711E+01 

0.4422E+01 

0. lOOOE+Ol 

699 

0 . 5726E+07 

8.155 

0.6910E+00 

0.6638E+00 

0 . 2974E+01 

0.4949E+01 

0.10008*01 

700 

0 . 5732E+07 

5.982 

0 . 6920E+00 

Q.6646E+Q0 

0.2981E+01 

0.4963E+01 

0.1000E*01 

725 

0.59888+07 

0.487 

0.7350E+00 

0.6966E+00 

0 . 3296E+0 1 

0 . 5592E+0 1 

0.10008*01 

745 

0 . 6254E+07 

3.103 

0.78228+00 

.0 .7300E+00 

0.3703E+01 

0.6407E+01 

0.10008+01 

Total  Transfer  Delta-V:  5.8264  ka /s  (  5826.39  a/s) 

Total  Transfer  Tiae:  72.394  days  (  0.20  yr) 
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